# Adaptive Neuro-Fuzzy Control of a Spherical Rolling Robot Using Sliding Mode
Control Theory-Based Online Learning Algorithm

Erkan Kayacan , Erdal Kayacan , Herman Ramon and Wouter Saeys E. Kayacan, E.
Kayacan, H. Ramon ans W. Saeys are with the Division of Mechatronics,
Biostatistics and Sensors, Department of Biosystems, KU Leuven, Kasteelpark
Arenberg 30, B-3001 Leuven, Belgium. e-mail: {erkan.kayacan, erdal.kayacan,
herman.ramon<EMAIL_ADDRESS>
###### Abstract

As a model is only an abstraction of the real system, unmodeled dynamics,
parameter variations and disturbances can result in poor performance of a
conventional controller based on this model. In such cases, a conventional
controller cannot remain well-tuned. This paper presents the control of a
spherical rolling robot by using an adaptive neuro-fuzzy controller in
combination with a SMC theory-based learning algorithm. The proposed control
structure consists of a neuro-fuzzy network and a conventional controller
which is used to guarantee the asymptotic stability of the system in a compact
space. The parameter updating rules of the fuzzy-neuro system using SMC theory
are derived, and the stability of the learning is proven using a Lyapunov
function. The simulation results show that the control scheme with the
proposed SMC theory-based learning algorithm is able to not only eliminate the
steady state error but also to improve the transient response performance of
the spherical rolling robot without knowing its dynamic equations.

##  I Introduction

Mobil robots are increasingly used in a variety of non-industrial applications
such as security surveillance, search and rescue, children education,
entertainment, etc. Spherical rolling mechanisms exhibit a number of
advantages with respect to wheeled and legged mechanisms. All mechanical and
electrical components including the actuation mechanism are securely located
inside a spherical shell rolling itself over the ground surface. The motion of
a sphere rolling without slipping over a surface is governed by non-holonomic
constraints, therefore spherical rolling robots are classified as non-
holonomic mechanical systems.

The fundamental difference between spherical and wheeled rolling motions is
the instantaneous number of degrees of freedom (DOF) between the mobile body
and ground surface. Since the sphere can simultaneously rotate around the
transverse and longitudinal axes, its instantaneous mobility is greater than
that of a wheel. Besides, spherical rolling mechanisms can change their
direction of motion easier than wheeled mechanisms. Unlike wheeled and legged
mechanisms, spheres can not fall over. The general problem of stability of
equilibrium frequently encountered in mobile robotics is naturally avoided
with the use of spherical rolling. On the other hand, highly complex nonlinear
equations are needed to describe the dynamics of spherical rolling systems.

A number of spherical robot prototypes have been developed in recent years.
One of the present mechanical structures consists of a wheeled vehicle located
at the bottom of the sphere. In such a scheme, the motion is provided by the
interaction of the vehicle’s wheels and the sphere. The wheeled vehicle’s
motion inside the sphere being also non-holonomic, the overall system is the
combination of two non-holonomic mechanical systems [1]. A linear model for
the longitudinal dynamics of the vehicle has been proposed to simplify the
equations of motion of the system [2]. Since spherical rolling robots have
highly nonlinear and very complex dynamics, the linearization approach is not
proper for the study of these systems. A similar design utilizes a single
wheel resting on the bottom of the sphere [3].

In another mechanical design, the driving motion of the rolling robot is
obtained by changing the position of the mass center inside the sphere, which
provides a gravitational torque [4, 5]. Such as design has been proposed in
[6, 7] where two DC motors associated with two flywheels are mounted inside
the sphere. When one of the flywheels turns around an arbitrarily chosen
reference direction, the spherical rolling robot rolls in the opposite
direction due to the conservation of angular momentum. However, such a
mechanism cannot roll over curvilinear trajectories. In order to change its
direction of motion, the sphere must stop rolling and then turn left or right
around the vertical axis. A similar design utilizes two perpendicular rotors
mounted inside the sphere [8].

A spherical robot with one DOF pendulum as actuation mechanism was designed in
[9]. A simplified dynamic model of a spherical robot with one DOF pendulum is
derived under some assumptions in [10, 11]. One of these assumptions is that
the longitudinal and lateral motions are decoupled from each other. Motion
control has been established by feedback linearization considering the
decoupled dynamics [12] and only driving motion [13].

Having different mechanical designs, spherical robots are used for very
different purposes. For example, Michaud [14, 15] invented a spherical mobile
robot called Roball, which is a tool in child-development studies. Roball, an
autonomous robot, contributes to the development of children’s language,
affective, intellectual and social skills. Longitudinal and lateral models of
Roball are discussed in [16].

Due to the highly nonlinear equations constituting the dynamic, a kinematic
model of the spherical rolling robots is usually used for path planning and in
the design of controllers. A trajectory was simulated using a kinematic model
in [17]. Under some assumptions, a feedback controller for a kinematic model
and a dynamic model based on the backstepping method are also designed [18].

In most of the papers in literature, some assumptions are made to obtain
simplified dynamic models for spherical rolling robots. For example, feedback
control and feedback linearization have been applied to the system by using
these simplified dynamic models [19]. Furthermore, the linearized dynamic
model is also obtained and some linear control methods are applied to the
system. However, these approaches are not feasible in practice due to the
highly complex nonlinear equations. Moreover, although the kinematic model has
been used to control the system because of its simplicity, it does not consist
of the system’s properties such as, masses, moment of inertia, etc… Thus,
kinematic models are not sufficient to control the spherical rolling robots.

In practice, the system behavior is affected by unmodeled dynamics, parameter
variations, uncertainties and disturbances. For instance, when the spherical
rolling robot rolls on a surface, the static and the dynamic friction cannot
be neglected. Since conventional controllers are time-invariant controllers,
these terms cause discontinuities and nonlinearities which render conventional
control invalid. Therefore, advanced intelligent control techniques are needed
to overcome the stated shortcomings. The control scheme used in this paper is
called feedback-error-learning and was firstly proposed in [20] for the
control of robot manipulators. It is based on the parallel work of a fuzzy
neural network (FNN)-based controller and a conventional feedback controller.

All systems are nonlinear and continuous in nature but, the data flow between
receiver and transmitter occurs in discrete-time domain in practice resulting
in some modeling errors. In [21], the modeling errors are estimated by using
an artificial neural network (ANN) estimator, and tracking errors are reduced
in discrete nonlinear systems. Even if the dynamic equations for the system at
hand would be available, uncertainties can exist because of the noise on the
sensors, environmental changes and/or nonlinear characteristics of the
actuators. Radial basis function neural networks with an online learning
algorithm have been proposed to handle the uncertainty problem [22]. The
dynamics of mobile robots with uncertainties and external disturbances have
been studied in [23]. Self-recurrent wavelet neural networks are also proposed
to estimate model uncertainties and external disturbances in the dynamics of
nonholonomic wheeled mobile robots [24]. Dual neuro-adaptive control for the
discrete-time dynamic control of non-holonomic mobile robots is introduced in
[25] which results in a major improvement in tracking performance, despite the
plant uncertainty and unmodeled dynamics.

As a solution to the modeling errors and uncertainties problems mentioned
above, an adaptive FNN-based controller with sliding mode control (SMC)
theory-based learning algorithm is proposed for the control of spherical
rolling robots. FNNs are a fusion of the capability of fuzzy reasoning to
handle uncertain information and the capability of ANNs to learn from input-
output data sets. Thus, FNNs are preferable approaches in engineering fields
[26, 27, 28]. Various learning algorithms have been proposed for FNNs. The
gradient-based algorithm works well when the system at hand has very slow
variations in its dynamics. However, since the gradient-based algorithms (e.g.
dynamic back propagation) include partial derivatives, the convergence speed
may be slow. Especially when the search space is complex. Moreover, the tuning
process can easily be trapped into a local minimum [29]. To alleviate the
problems mentioned, the use of evolutionary approaches has been suggested
[30]. However, the stability of such approaches is questionable and the
optimal values for the stochastic operators are difficult to derive.
Furthermore, the computational burden is very high. To overcome these issues,
a novel SMC theory-based algorithm is proposed in this study. The system can
provide adaptation to parameter variations, uncertainties and disturbances
through the SMC theory-based online learning algorithm and also shows robust
behaviour through the nature of SMC theory. SMC theory-based learning
algorithms cannot only make the overall system more robust, but also ensure
faster convergence than the traditional learning techniques in online tuning
of ANNs and FNNs [31]. The main contributions of this study beyond the state
of the art are a novel SMC theory-based online learning algorithm for FNNs and
the use of the proposed control algorithm in the velocity control of spherical
rolling robots. To the best knowledge of the authors, this is the first time
such an approach is ever used for the velocity control of a spherical rolling
robot.

This paper is organized as follows: The mathematical model of the spherical
rolling robot is presented in Section II. In Section III, the adaptive neuro-
fuzzy control approach is used to design an intelligent controller, and the
basics of the SMC theory-based online learning algorithm are given. Simulation
results and comparisons are presented in Section IV. Finally, a brief
conclusion of the study is given in Section V.

##  II The Mathematical Model of a Spherical Rolling Robot

The mathematical model of the spherical rolling robot is similar to the one
presented by Kayacan et. al. [19]. The main difference is that viscous
friction has been added into the equations of motion. A schematic illustration
of the spherical rolling robot is given in Fig. 1.

Figure 1: Modeling of rolling motion about transversal axis for overall
translation along $O-y$

###  II-A Kinematic Model

Formulation of the kinematic and dynamic equations governing the motion of the
spherical rolling mechanism with a pendulum is based on the following
assumptions:

  1. 1.

The sphere rolls over a perfectly horizontal surface without slipping.

  2. 2.

The center of mass of the entire system is also the geometric center of the
spherical shell.

  3. 3.

The pendulum is in vertical downward position when the sphere is in static
equilibrium.

A schematic illustration of the spherical rolling robot is given in Fig. 1.
$R_{f_{0}}$, $R_{f_{1}}$ and $R_{f_{2}}$ denote respectively, the reference
frames $O_{0}-X_{0}Y_{0}Z_{0}$, $O_{1}-X_{1}Y_{1}Z_{1}$ and
$O_{2}-X_{2}Y_{2}Z_{2}$. $R_{f_{0}}$ represents the inertial reference frame
fixed to the motion ground. $R_{f_{1}}$ is a moving frame attached to the
center of the sphere and allowed to translate only with respect to
$R_{f_{0}}$. $R_{f_{2}}$ is another moving frame attached also to the center
of the sphere but allowed to rotate only with respect to $R_{f_{1}}$ .
Relative angular positions between these frames can be described by several
methods such as Euler angles, Tait-Bryan angles, Roll-Pitch-Yaw angles, etc.

The variables of the spherical rolling robot are represented in Table I.

TABLE I: NOMENCLATURE $\theta$ | Rolling angle of the sphere around the x axis
---|---
$\alpha$ | Degree of freedom of pendulum
R | Radius of the sphere
l | Distance between the center of the sphere and
| the center of the pendulum
g | Gravitational acceleration
$\zeta$ | Damping coefficient

The angular velocity vector $\mathbf{\omega_{s}}$ and linear velocity vector
$\mathbf{v_{s}}$ of the center of the sphere, the position vector
$\mathbf{r_{p_{0}}}$, the angular velocity vector $\mathbf{\omega_{p}}$ and
the linear velocity vector $\mathbf{v}_{p}$ of the mass center of the pendulum
are given by:

$\displaystyle\mathbf{\omega_{s}}$ $\displaystyle=$
$\displaystyle-\dot{\theta}\mathbf{i}$ (1) $\displaystyle\mathbf{v_{s}}$
$\displaystyle=$ $\displaystyle-R\dot{\theta}\mathbf{j}$ (2)
$\displaystyle\mathbf{r_{p}}$ $\displaystyle=$ $\displaystyle
l\sin{(\alpha-\theta)}\mathbf{j}-l\cos{(\alpha-\theta)}\mathbf{k}$ (3)
$\displaystyle\mathbf{\omega_{p}}$ $\displaystyle=$
$\displaystyle(\dot{\alpha}-\dot{\theta})\mathbf{i}$ (4)
$\displaystyle\mathbf{v}_{p}$ $\displaystyle=$
$\displaystyle\big{(}-R\dot{\theta}+(\dot{\alpha}-\dot{\theta})l\cos{(\alpha-\theta)}\big{)}\mathbf{j}$
(5)
$\displaystyle+\big{(}(\dot{\alpha}-\dot{\theta})l\sin{(\alpha-\theta)}\big{)}\mathbf{k}$

where $\theta$, $\alpha$, $R$ and $l$ represent the rolling angle of the
sphere around the x axis, the rotation of the pendulum around the x axis, the
radius of the sphere and the distance between the center of the sphere and the
center of the pendulum, respectively. Similarly, $\mathbf{i}$, $\mathbf{j}$
and $\mathbf{k}$ represent the unit vector on the x, y and z axes,
respectively. It is to be noted that x,y and z axes denote the inertial
reference frame fixed to the motion ground.

###  II-B Dynamic Model

Let $E_{k}$ and $E_{p}$ denote respectively the total kinetic and potential
energy of the system. With $M_{s}$ and $m_{p}$ representing the masses of the
sphere and pendulum, $I_{s}$ and $I_{p}$ representing the mass moment of
inertia, and $v_{s}$, $\omega_{s}$, $v_{p}$ and $\omega_{p}$ representing the
linear and angular velocities of the sphere and pendulum, $r_{p-z}$
representing the vertical position of the mass center of the pendulum, the
Lagrangian function $L$ including only the terms due to rotations around the
transversal axis is then written as follows:

$\displaystyle L$ $\displaystyle=$ $\displaystyle E_{k}-E_{p}$ (6)
$\displaystyle=$
$\displaystyle\frac{1}{2}M_{s}\|\mathbf{v_{s}}\|^{2}+\frac{1}{2}I_{s}\|\mathbf{\omega_{s}}\|^{2}+\frac{1}{2}m_{p}\|\mathbf{v_{p}}\|^{2}+\frac{1}{2}I_{p}\|\mathbf{\omega_{p}}\|^{2}$
$\displaystyle-m_{p}gr_{p-z}$ $\displaystyle=$
$\displaystyle\frac{1}{2}M_{s}(-R\dot{\theta})^{2}+\frac{1}{2}I_{s}(-\dot{\theta})^{2}+\frac{1}{2}I_{p}(\dot{\alpha}-\dot{\theta})^{2}$
$\displaystyle+\frac{1}{2}m_{p}\Big{(}\big{(}-R\dot{\theta}+(\dot{\alpha}-\dot{\theta})l\cos{(\alpha-\theta)}\big{)}^{2}$
$\displaystyle+\big{(}(\dot{\alpha}-\dot{\theta})l\sin{(\alpha-\theta)}\big{)}^{2}\Big{)}-m_{p}gl\cos{(\alpha-\theta)}$

It is assumed that the viscous friction operates between the sphere and the
surface. The loss due to the viscous friction is written in an energy
dissipation function that depends on the velocities of the system and the
damping constant:

$S=\frac{1}{2}\zeta\dot{q}_{i}^{2}=\frac{1}{2}\zeta(\dot{\theta}^{2}+\dot{\alpha}^{2})$
(7)

For translation along $O-y$, the Euler-Lagrange equations of the system are
written as follows:

$\frac{d}{dt}(\frac{\partial L}{\partial\dot{q}_{i}})-\frac{\partial
L}{\partial q_{i}}+\frac{\partial S}{\partial\dot{q}_{i}}=Q_{i}$ (8)

where $q_{1}=\theta$ and $q_{2}=\alpha$ are the generalized coordinates. In
fact, when the pendulum is rotated through an input torque, a reaction torque
about the shaft occurs in the opposite direction [9]. $Q_{1}=Q_{2}=\tau$
represents the input torque to rotate the pendulum.

$\displaystyle\begin{array}[]{c}Q_{\theta}=\tau\\\ Q_{\alpha}=\tau\\\
\end{array}$ (11)

As the equations of motion of a mechanical system can be written as follows:

$M\big{(}q(t)\big{)}\ddot{q(t)}+C\big{(}q(t),\dot{q(t)}\big{)}+G\big{(}q(t)\big{)}=u(t)$
(12)

The equations of motion can be finally written in the following matrix form:

$\displaystyle\left[\begin{array}[]{cc}M_{11}&M_{12}\\\ M_{21}&M_{22}\\\
\end{array}\right]\left[\begin{array}[]{c}\ddot{\theta}\\\ \ddot{\alpha}\\\
\end{array}\right]+\left[\begin{array}[]{c}C_{11}\\\ C_{21}\\\
\end{array}\right]+\left[\begin{array}[]{c}G_{11}\\\ G_{21}\\\
\end{array}\right]=\left[\begin{array}[]{c}\tau\\\ \tau\\\ \end{array}\right]$
(23)

where

$\displaystyle M_{11}$ $\displaystyle=$ $\displaystyle
M_{s}R^{2}+m_{p}R^{2}+m_{p}l^{2}+I_{s}+I_{p}$
$\displaystyle+2m_{p}Rl\cos{(\alpha-\theta)}$ $\displaystyle M_{12}$
$\displaystyle=$ $\displaystyle
M_{21}=-m_{p}l^{2}-I_{p}-m_{p}Rl\cos{(\alpha-\theta)}$ $\displaystyle M_{22}$
$\displaystyle=$ $\displaystyle m_{p}l^{2}+I_{p}$ $\displaystyle C_{11}$
$\displaystyle=$ $\displaystyle
m_{p}Rl\sin{(\alpha-\theta)}(\dot{\alpha}-\dot{\theta})^{2}+\zeta\dot{\theta}$
$\displaystyle C_{21}$ $\displaystyle=$ $\displaystyle\zeta\dot{\alpha}$
$\displaystyle G_{11}$ $\displaystyle=$ $\displaystyle-
G_{21}=-m_{p}gl\sin{(\alpha-\theta)}$

In this study, velocity control of the spherical rolling robot is studied by
using the above equations of motion. If rolling motion over curvilinear
trajectories is desired, the equations of motion presented in [19] can be
used. Different from the approach in [19], a viscous friction term is added in
this study.

##  III The Adaptive Neuro-Fuzzy Control Approach

###  III-A The Control Scheme and the Adaptive Neuro-Fuzzy Inference System

In real-time applications, structural and parametric uncertainties such as
unmodeled dynamics and physical disturbances cause unwanted effects on the
system behaviour. Like tire dynamics, the spherical rolling robot dynamics are
affected by the slip ratio and variations of the surface conditions. Thus,
spherical rolling robots should have an adaptive intelligent controller with
respect to the high uncertainty and the variation in the surface properties.
While conventional controllers struggle with these problems, neuro-fuzzy
control can overcome these limitations and provide higher robustness.

In the presented approach, a conventional PD or PID controller and a neuro-
fuzzy controller are working in parallel. The conventional controller is used
to guarantee the global asymptotic stability of the system in compact space.
The PD control law is written as follows:

$\tau_{c}=k_{p}e+k_{d}\dot{e}\\\ $ (24)

where $e=\theta_{d}-\theta$ is the feedback error, $\theta_{d}$ is the desired
position value, $k_{p}$ and $k_{d}$ are the controller gains.

The proposed control scheme for the PD controller case, the PID controller
case and the structure of the proposed FNN are illustrated in Fig. 2, Fig. 3
and Fig. 4, respectively. The terms of $\alpha$ and $\beta$ are the
coefficients of the PI controller. As can be seen from the proposed control
scheme, a FNN which has two inputs $x_{1}(t)=e(t)$, $x_{2}(t)=\dot{e}(t)$ and
one output $f_{ij}$, is applied as a feedback controller.

Figure 2: Block diagram of the proposed control scheme for the PD controller
case Figure 3: Block diagram of the proposed control scheme for the PID
controller case Figure 4: The FNN

The fuzzy _If-Then_ rule of a zeroth-order Takagi-Sugeno-Kang (TSK) model with
two input variables where the consequent part is a linear function of the
input variables can be defined as follows:

$R_{ij}:\;\;\textrm{If}\;x_{1}\;\textrm{is}\;\;A_{i}\;\;\textrm{and}\;x_{2}\;\textrm{is}\;\;B_{j},\;\;\textrm{Then}\;f_{ij}=d_{ij}$
(25)

where $x_{1}$ and $x_{2}$ are the inputs of the FNN, $A_{i}$ and $B_{j}$ are
fuzzy sets corresponding to the input 1 and input 2, respectively. The zeroth-
order function, $f_{ij}$ is the consequent part of the rules where $I$
$(i=1,\dots,I)$ and $J$ $(j=1,\dots,J)$ being the number of membership
functions used for the input 1 and input 2, respectively.

The following are the basic explanations for Fig. 4:

•Layer 1: The first layer is called input later. This layer maps the crisp
inputs $x_{1}$ and $x_{2}$ into fuzzified values using membership functions.
Every node has a node function described as follows:

$\displaystyle O^{1}_{i}=\mu_{A_{i}}(x_{1})$ $\displaystyle
O^{1}_{j}=\mu_{B_{j}}(x_{2})$ (26)

where $\mu_{A_{i}}(x_{1})$ and $\mu_{B_{j}}(x_{2})$ are the membership values
for the inputs 1 and 2, respectively. In the case of Gaussian membership
functions, $\mu_{A_{i}}(x_{1})$ and $\mu_{B_{j}}(x_{2})$ are written as
follows:

$\displaystyle\mu_{A_{i}}(x_{1})=\exp\Big{(}-(\frac{x_{1}-c_{A_{i}}}{\sigma_{A_{i}}})^{2}\Big{)}$
$\displaystyle\mu_{B_{j}}(x_{2})=\exp\Big{(}-(\frac{x_{2}-c_{B_{j}}}{\sigma_{B_{j}}})^{2}\Big{)}$
(27)

where $c$ and $\sigma$ are respectively the mean and the standard deviation of
the membership functions. These parameters $c,\sigma>0$ are the tunable
parameters of the neuro-fuzzy structure.

•Layer 2: The strength of the rule is obtained in this layer. The output of
each node, which calculates the firing strength $w_{ij}$ of a rule, is
multiplied with all incoming signals. The node function of this layer is
written as follows:

$O^{2}_{ij}=w_{ij}=\mu_{A_{i}}(x_{1})\mu_{B_{j}}(x_{2})$ (28)

The previous equation (28) is rewritten considering (III-A) as follows:

$O^{2}_{ij}=w_{ij}=\exp\Big{(}-(\frac{x_{1}-c_{A_{i}}}{\sigma_{A_{i}}})^{2}-(\frac{x_{2}-c_{B_{j}}}{\sigma_{B_{j}}})^{2}\Big{)}$
(29)

•Layer 3: The normalized firing strength of each node is written as follows:

$O^{3}_{ij}=\overline{w}_{ij}=\frac{w_{ij}}{\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}{w}_{ij}}$
(30)

•Layer 4: Each node in this layer is an adaptive node with node function
written by the following equation:

$O^{4}_{ij}=\overline{w}_{ij}f_{ij}$ (31)

•Layer 5: There is one node where all incoming signals are summed in this
layer. The output signal of the neuro-fuzzy network $\tau_{n}$ is calculated
by the function described as follows:

$O^{5}=\tau_{n}(t)=\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\overline{w}_{ij}f_{ij}=\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\overline{w}_{ij}d_{ij}$
(32)

The control input fed to the system is the overall torque $\tau$ determined as
follows:

$\tau=\tau_{c}-\tau_{n}$ (33)

where $\tau_{c}$ and $\tau_{n}$ are the input torques generated by the PD
controller and the neuro-fuzzy feedback controller, respectively.

###  III-B The Sliding Mode Learning Algorithm

The input signals and their time derivatives
$x_{1}(t),x_{2}(t),\dot{x}_{1}(t),\dot{x}_{2}(t)$ cannot have infinite values.
Thus, they are bounded as follows:

$\displaystyle\mid x_{1}(t)\mid\leq B_{x},\;\;\mid x_{2}(t)\mid\leq
B_{x}\quad\forall t$ $\displaystyle\mid\dot{x}_{1}(t)\mid\leq
B_{\dot{x}},\;\;\mid\dot{x}_{2}(t)\mid\leq B_{\dot{x}}\quad\forall t$

Similarly, the input torque to the system and its time derivative are bounded:

$\mid\tau(t)\mid\leq B_{\tau},\;\;\mid\dot{\tau}(t)\mid\leq
B_{\dot{\tau}}\quad\forall t$ (35)

The output of the PD controller $\tau_{c}\left(t\right)$ is defined as a time-
varying sliding surface by using the principles of SMC theory [32]:

$S_{c}\left(\tau_{n},\tau\right)=\tau_{c}\left(t\right)=\tau_{n}\left(t\right)+\tau\left(t\right)=0$
(36)

The sliding surface $S_{p}\left(e,\dot{e}\right)$ is defined as follows:

$S_{p}\left(e,\dot{e}\right)=\dot{e}+\lambda e$ (37)

where $\lambda$ is a constant determining the slope of the sliding surface.

Definition: A sliding motion exists on the sliding surface
$S_{c}\left(\tau_{n},\tau\right)=\tau_{c}\left(t\right)=0$ after time $t_{h}$,
if the condition
$S_{c}(t)\dot{S}_{c}(t)=\tau_{c}\left(t\right)\dot{\tau}_{c}\left(t\right)<0$
is satisfied for all $t$ in some nontrivial semi-open subinterval of time of
the form $\left[t,t_{h}\right)\subset\left(-\infty,t_{h}\right)$.

Since it is desired to design a dynamical feedback adaptation mechanism, or
online learning algorithm for the FNN parameters, the sliding mode condition
of the above definition is applied.

######  Theorem 1.

The adaptation laws for the parameters of the proposed FNN are given by the
following equations:

$\dot{c}_{A_{i}}=\dot{x_{1}}+s_{A_{i}}\alpha\textrm{sgn}\left(\tau_{c}\right)$
(38)

$\dot{c}_{B_{j}}=\dot{x_{2}}+s_{B_{j}}\alpha\textrm{sgn}\left(\tau_{c}\right)$
(39)

$\dot{\sigma}_{A_{i}}=-\bigg{(}\sigma_{A_{i}}+\frac{(\sigma_{A_{i}})^{3}}{s_{A}^{T}s_{A}}\bigg{)}\alpha\textrm{sgn}\left(\tau_{c}\right)$
(40)

$\dot{\sigma}_{B_{j}}=-\bigg{(}\sigma_{B_{j}}+\frac{(\sigma_{B_{j}})^{3}}{s_{B}^{T}s_{B}}\bigg{)}\alpha\textrm{sgn}\left(\tau_{c}\right)$
(41)

$\dot{f}_{ij}=-\frac{\overline{W}_{ij}}{\overline{W}^{T}\overline{W}}\alpha\textrm{sgn}\left(\tau_{c}\right)$
(42)

The vectors in (38)-(41) are defined as follows:

$\displaystyle
s_{A_{i}}=x_{1}-c_{A_{i}}\;\textrm{and}\;s_{A}=\left[s_{A_{1}}\,s_{A_{2}}...\,s_{A_{I}}\right]^{T}$
(43) $\displaystyle
s_{B_{j}}=x_{2}-c_{B_{j}}\;\textrm{and}\;s_{B}=\left[s_{B_{1}}\,s_{B_{2}}...\,s_{B_{J}}\right]^{T}$
(44)

where $\alpha$ is the learning rate which is a sufficiently large positive
design constant which satisfies the following inequality:

$B_{\dot{\tau}}<\alpha$ (45)

This ensures that for a given arbitrary initial condition $\tau_{c}(0)$, the
learning error $\tau_{c}(t)$ will converge to zero within finite time $t_{h}$.

_Proof:_ The reader is referred to Appendix A.

If $\lambda$ is considered $\lambda=\frac{k_{p}}{k_{d}}$, the relation between
the sliding line $S_{p}$ and the zero adaptive learning error level $S_{c}$
can be written as follows:

$S_{c}=\tau_{c}=k_{d}\dot{e}+k_{p}e=k_{d}\left(\dot{e}+\frac{k_{p}}{k_{d}}e\right)=k_{d}S_{p}$
(46)

The tracking performance of the velocity control system in the spherical
rolling robot is analyzed by the following Lyapunov function candidate:

$V_{p}=\frac{1}{2}S_{p}^{2}$ (47)

######  Theorem 2.

If the adaptation strategy for the adjustable parameters of the FNN is chosen
as in (38)-(42), then the negative definiteness of the time derivative of the
Lyapunov function in (47) is ensured.

_Proof:_ The reader is referred to Appendix B.

##  IV Simulation Results and Discussion

The numerical values used in [4, 6, 9] are considered for the numerical values
in this study which are $M_{s}=3$ $kg$, $m_{p}=2$ $kg$, $R=0.2$ $m$, $l=0.075$
$m$ and $g=9.81$ $m/s^{2}$. The damping coefficient, $\zeta$, in the equations
of motion is set to $0.2$, and the sampling period of the simulations is set
to $0.001$ $s$. The number of membership functions for the input $1$ and input
$2$ is chosen as $I=J=3$ for all the simulations.

SMC theory suffers from high-frequency oscillations called _chattering_
because of the input control law. Several approaches have been suggested to
get rid of this problem. In this paper, the $sgn$ function in (38)-(42) is
replaced by the following equation to decrease the chattering effect:

$\textrm{sgn}\left(\tau_{c}\right):=\frac{\tau_{c}}{\mid\tau_{c}\mid+\delta}$
(48)

where $\delta=0.05$.

###  IV-A Case 1: PD controller

The coefficients of the PD controller are set to $k_{p}=1$ and $k_{d}=0.05$ by
trial-and-error method. The following reference signal has been applied to the
system (49):

$\textrm{Reference}(t)=\left\\{\begin{array}[]{l
l}1\;\textrm{rad/s}&\quad\textrm{if $0<t\leq 5$ }\\\
2\;\textrm{rad/s}&\quad\textrm{if $5<t\leq 10$ }\\\
1.5\;\textrm{rad/s}&\quad\textrm{if $10<t\leq 15$ }\\\ \end{array}\right.$
(49)

Figures 5-6 show the velocity response and the error of the system for the PD
controller working alone and the PD controller working in parallel with the
FNN, respectively. As can be seen from Fig. 5, the PD controller cannot
eliminate the steady state error for a time-varying step input. The system
does not have any steady state error for the case of the FNN working in
parallel with a PD controller.

Figure 5: The velocity response of the system for PD and PD+FNN controllers
Figure 6: The error of the system for PD and PD+FNN controllers

###  IV-B Case 2: PID controller

The coefficients of the PID controller are set to
$k_{p}=1,k_{d}=0.05,\alpha=1$ and $\beta=2$ by trial-and-error method. The
reference signal for this simulation is the same the one (See 49) in Case 1.

In Figure 7 the responses of the PID controller and the combination of the FNN
in parallel with this PID controller are illustrated. The FNN learns the
system dynamics after a finite time duration, and results in a smaller rise
time, overshoot and settling time. From Fig. 7 and 8 it can be concluded that
the control scheme with the proposed SMC theory-based learning algorithm is
able to not only eliminate the steady state error (shown in Case 1) but also
to improve the transient response performance of the system (shown in Case 2).
Figure 8 shows the zoomed view of Fig. 7 between $0^{th}-3^{rd}$,
$5^{th}-8^{th}$ and $10^{th}-13^{th}$.

Figure 7: The velocity response of the system for PID and PID+FNN controllers
Figure 8: Zoomed view of Figure 7

To test the robustness of this approach, the damping coefficient is set to
$0.5$ with a noise level $SNR=20dB$. In Figure 9 the responses of the PID
controller and the combination of the FNN working in parallel with the PID
controller are illustrated. In Fig. 10 close ups of Fig. 9 are shown. It can
be observed that the adaptive neuro-fuzzy control approach is much more robust
to the uncertainties and gives more satisfactory results regarding smaller
rise time and settling time than the PID controller working alone.

Figure 9: The velocity response of the system for PID and PID+FNN controllers
when the damping coefficient is set to $0.5$ with a noise level $SNR=20dB$
Figure 10: Zoomed view of Figure 9

In a next step, the following variation of the damping coefficient over time
was applied to the system to test the robustness of the controllers:

$\zeta(t)=\left\\{\begin{array}[]{l l}0.2\;\textrm{rad/s}&\quad\textrm{if
$0<t\leq 5$ }\\\ 0.5\;\textrm{rad/s}&\quad\textrm{if $5<t\leq 10$ }\\\
0.8\;\textrm{rad/s}&\quad\textrm{if $10<t\leq 15$ }\\\ \end{array}\right.$
(50)

In Figure 11 the velocity responses of a conventional PID controller and the
combination of the FNN working in parallel with this PID controller are shown,
while close ups are shown in Fig. 12. As can be seen from Fig. 11 \- 12, the
FNN can adapt its parameters when the coefficient of viscous friction changes
suddenly. For the case when the coefficient of the viscous friction changes,
while the settling time for the PID controller working alone is approximately
$2.5$ seconds, it is less than $1$ second for the FNN working in parallel with
a PID controller. From these simulations it can be concluded that the adaptive
neuro-fuzzy control approach is much more robust to the parameter variations
compared to the case of the PID controller working alone.

Figure 11: The velocity response of the system for PID and PID+FNN controllers
when the damping coefficient is set to $0.2$, $0.5$ and $0.8$ at $0^{th}$
second, $5^{th}$ second and $10^{th}$ second, respectively Figure 12: Zoomed
view of 11

In Figure 13 the control signals coming from the PID controller alone and the
FNN working in parallel with the PID controller are illustrated for the case
when the damping coefficient is set to $0.2$, $0.5$ and $0.8$ at $0^{th}$
second, $5^{th}$ second and $10^{th}$ second, respectively. In order not to
confuse the reader, since the control input to the system is equal to
$\tau_{c}-\tau_{n}$ when the PID controller works in parallel with the FNN,
$\tau_{n}$ is multiplied by $-1$ in this figure. As can be seen from Fig. 13,
the output of the PID controller is approximately zero while the output of the
FNN changes as the damping coefficient changes as well. In other words, the
FNN takes the responsibility for controlling the system after a finite time
duration.

Figure 13: The control signals coming from the PID controller and the FNN when
the damping coefficient is set to $0.2$, $0.5$ and $0.8$ at $0^{th}$ second,
$5^{th}$ second and $10^{th}$ second, respectively

##  V Conclusion

In this study, an adaptive neuro-fuzzy controller with SMC theory-based online
learning has been elaborated for velocity control of a spherical rolling robot
subject to parameter variations and uncertainties in its dynamics. The
simulation studies show that the proposed adaptive neuro-fuzzy control scheme
(a conventional controller working in parallel with the FNN) results in a
better performance and higher robustness than when compared to the
conventional stand-alone controller. The proposed control algorithm is able to
not only eliminate the steady state error in the case of conventional PD
stand-alone controller but also improve the transient response performance of
the system in the case of conventional PID stand-alone controller. Thanks to
the SMC theory-based online learning algorithm the parameters of the
controller are automatically adapted to cope with the parameter variations and
uncertainties. In addition to its robustness, the control approach with the
proposed learning algorithm is computationally simple, especially when
compared to the gradient descent and the evolutionary algorithms.

##  Appendix A Proof of Theorem 1

The time derivative of (III-A) is written as follows:

$\displaystyle\dot{\mu}_{A_{i}}(x_{1})=-2N_{A_{i}}(N_{A_{i}})^{{}^{\prime}}\mu_{A_{i}}(x_{1})$
$\displaystyle\dot{\mu}_{B_{j}}(x_{2})=-2N_{B_{j}}(N_{B_{j}})^{{}^{\prime}}\mu_{B_{j}}(x_{2})$
(51)

where

$\displaystyle N_{A_{i}}=\frac{x_{1}-c_{A_{i}}}{\sigma_{A_{i}}}$
$\displaystyle N_{B_{j}}=\frac{x_{2}-c_{B_{j}}}{\sigma_{B_{j}}}$ (52)

The time derivative of (A) is written as follows:

$\displaystyle\dot{N}_{A_{i}}=\frac{(\dot{x}_{1}-\dot{c}_{A_{i}})\sigma_{A_{i}}-(x_{1}-c_{A_{i}})\dot{\sigma}_{A_{i}}}{\sigma_{A_{i}}^{2}}$
$\displaystyle\dot{N}_{B_{j}}=\frac{(\dot{x}_{2}-\dot{c}_{B_{j}})\sigma_{B_{i}}-(x_{2}-c_{B_{j}})\dot{\sigma}_{B_{j}}}{\sigma_{B_{j}}^{2}}$
(53)

Combination of (A) and (A) gives:

$N_{A_{i}}\dot{N}_{A_{i}}=N_{B_{j}}\dot{N}_{B_{j}}=\alpha\textrm{sgn}\left(\tau_{c}\right)$
(54)

The time derivative of (30) is written as follows:

$\dot{\overline{w}}_{ij}=-\overline{w}_{ij}\dot{k}_{ij}+\overline{w}_{ij}\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\overline{w}_{ij}\dot{k}_{ij}$
(55)

where

$\dot{k}_{ij}=2\Big{(}N_{A_{i}}(N_{A_{i}})^{{}^{\prime}}+N_{B_{j}}(N_{B_{j}})^{{}^{\prime}}\Big{)}$
(56)

The stability of the proposed control approach is investigated by using the
following Lyapunov function:

$V=\frac{1}{2}\left(\tau_{c}\right)^{2}$ (57)

The time derivative of (57) is written as follows:

$\dot{V}=\dot{\tau}_{c}\tau_{c}=\tau_{c}(\dot{\tau}_{n}+\dot{\tau})$ (58)

where

$\dot{\tau}_{n}=\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}(\dot{f}_{ij}\overline{w}_{ij}+f_{ij}\dot{\overline{w}}_{ij})$
(59)

Substitution of (59), (55), (56), (54) and (42) into (58), gives:

$\displaystyle\dot{V}$ $\displaystyle=$
$\displaystyle\tau_{c}\Big{[}\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\Big{(}\dot{f}_{ij}\overline{w}_{ij}+f_{ij}(-\overline{w}_{ij}\dot{k}_{ij}$
(60) $\displaystyle+$
$\displaystyle\overline{w}_{ij}\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\overline{w}_{ij}\dot{k}_{ij})\Big{)}+\dot{\tau}\Big{]}$
$\displaystyle=$
$\displaystyle\tau_{c}\Bigg{[}\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\Bigg{(}\dot{f}_{ij}\overline{w}_{ij}$
$\displaystyle+$ $\displaystyle
f_{ij}\Big{(}-2\overline{w}_{ij}(N_{A_{i}}(N_{A_{i}})^{{}^{\prime}}+N_{B_{j}}(N_{B_{j}})^{{}^{\prime}})$
$\displaystyle+$ $\displaystyle
2\overline{w}_{ij}\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\overline{w}_{ij}(N_{A_{i}}(N_{A_{i}})^{{}^{\prime}}+N_{B_{j}}(N_{B_{j}})^{{}^{\prime}})\Big{)}\Bigg{)}+\dot{\tau}\Bigg{]}$
$\displaystyle=$
$\displaystyle\tau_{c}\Bigg{[}\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\Bigg{(}\dot{f}_{ij}\overline{w}_{ij}+f_{ij}\Big{(}-4\overline{w}_{ij}\alpha\textrm{sgn}\left(\tau_{c}\right)$
$\displaystyle+$ $\displaystyle
4\overline{w}_{ij}\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\overline{w}_{ij}\alpha\textrm{sgn}\left(\tau_{c}\right)\Big{)}\Bigg{)}+\dot{\tau}\Bigg{]}$
$\displaystyle=$
$\displaystyle\tau_{c}\Big{[}\sum\limits_{i=1}^{I}\sum\limits_{j=1}^{J}\dot{f}_{ij}\overline{w}_{ij}+\dot{\tau}\Big{]}$
$\displaystyle=$
$\displaystyle\tau_{c}\Big{[}-\alpha\textrm{sgn}\left(\tau_{c}\right)+\dot{\tau}\Big{]}$
$\displaystyle=$
$\displaystyle\Big{[}-\alpha\mid\tau_{c}\mid+\dot{\tau}\mid\tau_{c}\mid\Big{]}$

The time derivative $\dot{V}$ of the Lyapunov function $V$ must be smaller
than zero to satisfy the stability of the learning.

$\dot{V}=\Big{[}-\alpha\mid\tau_{c}\mid+\dot{\tau}\mid\tau_{c}\mid\Big{]}<0$
(61)

When $\dot{\tau}$ reaches its maximal value $B_{\dot{\tau}}$, (61) can be re-
written as follows:

$\dot{V}=\Big{[}-\alpha\mid\tau_{c}\mid+B_{\dot{\tau}}\mid\tau_{c}\mid\Big{]}<0\;\;\;\;if\;\;\;\;B_{\dot{\tau}}<\alpha$
(62)

##  Appendix B Proof of Theorem 2

The time derivative of the Lyapunov function in (47) is written as follows:

$\displaystyle\dot{V}_{p}$ $\displaystyle=$
$\displaystyle\dot{S}_{p}S_{p}=\frac{1}{k^{2}_{d}}\dot{S}_{c}S_{c}=\frac{1}{k^{2}_{d}}\dot{\tau}_{c}\tau_{c}$
(63) $\displaystyle=$ $\displaystyle\frac{1}{k^{2}_{d}}\dot{V}<0,\;\;\forall
S_{c},S_{p}\neq 0$

## References

  * [1] J. Alves and J. Dias, “Design and control of a spherical mobile robot,” _Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering_ , vol. 217, pp. 457–467, 2003.
  * [2] C. Camicia, F. Conticelli, and A. Bicchi, “Nonholonomic kinematics and dynamics of the sphericle,” in _Proceedings of IEEE/RSJ International Conference on Intelligent Robots and Systems_ , vol. 1, Takamatsu, Japan, 2000, pp. 805 – 810.
  * [3] A. Halme, T. Schonberg, and Y. Wang, “Motion control of a spherical mobile robot,” in _Proceedings of fourth International Workshop on Advanced Motion Control_ , vol. 1, Albuquerque, NM, USA, 1996, pp. 259–264.
  * [4] A. A. H. Javadi and P. Mojabi, “Introducing glory: A novel strategy for an omnidirectional spherical rolling robot,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 126, pp. 678–683, 2004.
  * [5] R. Mukherjee, M. A. Minor, and J. T. Pukrushpan, “Motion planning for a spherical mobile robot: Revisiting the classical ball-plate problem,” _ASME Journal of Dynamic Systems Measurement and Control_ , vol. 124, pp. 502–511, 2002.
  * [6] V. A. Joshi and R. N. Banavar, “Motion analysis of a spherical mobile robot,” _Robotica_ , vol. 27, pp. 343 – 353, 2009.
  * [7] V. A. Joshi, R. N. Banavar, and R. Hippalgaonkar, “Design and analysis of a spherical mobile robot,” _Mechanism and Machine Theory_ , 2009.
  * [8] S. Bhattacharya and S. K. Agrawal, “Spherical rolling robot: a design and motion planning studies,” _IEEE Transactions on Robotics and Automation_ , vol. 16, pp. 835–839, 2000.
  * [9] Y. Ming, D. Zongquan, Y. Xinyi, and Y. Weizhen, “Introducing HIT spherical robot: Dynamic modeling and analysis based on decoupled subsystem,” in _Proceedings of the IEEE International Conference on Robotics and Biomimetics_ , Harbin, China, 2006, pp. 181–186.
  * [10] Q. Jia, H. Sun, and D. Liu, “Analysis of actuation for a spherical robot,” in _Proceedings of the 2008 IEEE Conference on Robotics, Automation and Mechatronics_ , Chengdu, 2008, pp. 266–271.
  * [11] D. Liu, H. Sun, and Q. Jia, “Stabilization and path following of a spherical robot,” in _Proceedings of the 2008 IEEE Conference on Robotics, Automation and Mechatronics_ , Chengdu, 2008, pp. 676 – 682.
  * [12] D. Liu, H. Sun, Q. Jia, and L. Wang, “Motion control of a spherical mobile robot by feedback linearization,” in _Proceedings of the 7th World Congress on Intelligent Control and Automation_ , Chongqing, China, 2008, pp. 965–970.
  * [13] D. Liu, H. Sun, and Q. Jia, “A family of spherical mobile robot: Driving ahead motion control by feedback linearization,” in _Proceedings of 2nd International Symposium on Systems and Control in Aerospace and Astronautics_ , Shenzhen, China, 2008.
  * [14] F. Michaud, J.-F. Laplante, H. Larouche, A. Duquette, S. Caron, D. Letourneau, and P. Masson, “Autonomous spherical mobile robot for child-development studies,” _IEEE Transactions on Systems, Man and Cybernetics, Part A: Systems and Humans_ , vol. 35, pp. 471 – 480, 2005.
  * [15] F. Michaud and S. Caron, “Roball, the rolling robot,” _Autonomous Robots_ , vol. 12, pp. 211 – 222, 2002.
  * [16] J.-F. Laplante, P. Masson, and F. Michaud, “Analytical longitudinal and lateral models of a spherical rolling robot,” Laborius, Department of Electrical and Computer Engineering, Universite de Shebroke, Canada, Tech. Rep., 2002.
  * [17] H. Alizadeh and M. Mahjoob, “Effect of incremental driving motion on a vision-based path planning of a spherical robot,” in _Proceedings of the 2009 Second International Conference on Computer and Electrical Engineering_ , vol. 1, Dubai, United Arab Emirates, 2009, pp. 299 – 303.
  * [18] T. Otani, T. Urakuboi, S. Maekawa, H. Tamaki, and Y. Tada, “Position and attitude control of a spherical rolling robot equipped with a gyro,” in _Proceedings of the 9th IEEE International Workshop on Advanced Motion Control_ , Istanbul, 2006, pp. 416 – 421.
  * [19] E. Kayacan, Z. Y. Bayraktaroglu, and W. Saeys, “Modeling and control of a spherical rolling robot: A decoupled dynamics approach,” _Robotica_.
  * [20] H. Gomi and M. Kawato, “Neural network control for a closed-loop system using feedback-error-learning,” _Neural Networks_ , vol. 6, no. 7, pp. 933 – 946, 1993. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S089360800980004X
  * [21] D. Munoz and D. Sbarbaro, “An adaptive sliding-mode controller for discrete nonlinear systems,” _IEEE Transactions on Industrial Electronics_ , vol. 47, pp. 574 – 581, 2000.
  * [22] S. J. Huang, K. S. Huang, and K. C. Chiou, “Development and application of a novel radial basis function sliding mode controller,” _Mechatronic_ , vol. 13, pp. 313–329, 2003.
  * [23] T. Das and I. N. Kar, “Design and implementation of an adaptive fuzzy logic-based controller for wheeled mobile robots,” _IEEE Transactions on Control Systems Technology_ , vol. 14, pp. 501 – 510, 2006.
  * [24] B. S. Park, S. J. Yoo, J. B. Park, and Y. H. Choi, “Adaptive neural sliding mode control of nonholonomic wheeled mobile robots with model uncertainty,” _IEEE Transactions on Control Systems Technology_ , vol. 17, pp. 207 – 214, 2009.
  * [25] M. Bugeja, S. Fabri, and L. Camilleri, “Dual adaptive dynamic control of mobile robots using neural networks,” _Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on_ , vol. 39, no. 1, pp. 129 –141, feb. 2009.
  * [26] A. V. Topalov, Y. Oniz, E. Kayacan, and O. Kaynak, “Neuro-fuzzy control of antilock braking system using sliding mode incremental learning algorithm,” _Neurocomputing_ , vol. 74, no. 11, pp. 1883 – 1893, 2011.
  * [27] E. Kayacan, Y. Oniz, A. C. Aras, O. Kaynak, and R. Abiyev, “A servo system control with time-varying and nonlinear load conditions using type-2 tsk fuzzy neural system,” _Applied Soft Computing_ , vol. In Press, Corrected Proof, pp. –, 2011. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1568494611001049
  * [28] M. Ahmadieh Khanesar, E. Kayacan, M. Teshnehlab, and O. Kaynak, “Extended kalman filter based learning algorithm for type-2 fuzzy logic systems and its experimental evaluation,” _Industrial Electronics, IEEE Transactions on_ , vol. PP, no. 99, p. 1, 2011.
  * [29] A. Topalov and O. Kaynak, “Online learning in adaptive neurocontrol schemes with a sliding mode algorithm,” _IEEE Transactions on Systems, Man, and Cybernetics – Part B: Cybernetics_ , vol. 31, no. 3, pp. 445 –450, december 2001\.
  * [30] Y. Zhou and M. Er, “An evolutionary approach toward dynamic self-generated fuzzy inference systems,” _Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on_ , vol. 38, no. 4, pp. 963 –969, aug. 2008.
  * [31] M. O. Efe, O. Kaynak, and X. Yu, “Sliding mode control of a three degrees of freedom anthropoid robot by driving the controller parameters to an equivalent regime,” _ASME Journal of Dynamic Systems, Measurement, and Control_ , vol. 122, no. 4, pp. 632 –640, December 2000.
  * [32] V. I. Utkin, _Sliding Modes in Control Optimization_. Springer-Verlag, 1992.

|  Erkan Kayacan (S’12) was born in Istanbul, Turkey, on April 17, 1985. He
received the B.Sc. and the M.Sc. degrees in mechanical engineering from
Istanbul Technical University, Istanbul, in 2008 and 2010, respectively. He is
currently working toward the Ph.D. degree in the Department of Biosystems, KU
Leuven. He is currently a Research Assistant in KU Leuven at the Department of
Biosystems (BIOSYST) in the Division of Mechatronics, Biostatistics and
Sensors (MeBioS). His research interests include robotics, nonlinear control,
intelligent control, system identification, fuzzy theory, and grey system
theory.
---|---
|  Erdal Kayacan (S’06-SM’12) was born in Istanbul, Turkey on January 7, 1980.
He received the B.Sc. degree in electrical engineering from Istanbul Technical
University, Istanbul, Turkey, in 2003. He received the M.Sc. and Ph.D. degrees
in systems and control engineering and electrical and electronics engineering
from Bogazici University, Istanbul, Turkey, in 2006 and 2011, respectively. He
is currently a post-doctoral researcher in KU Leuven at the department of
biosystems (BIOSYST) in the division of mechatronics, biostatistics and
sensors (MeBioS). His research interests include large scale systems, soft
computing, intelligent control, fuzzy logic theory. Dr. Kayacan is active in
IEEE CIS Student Activities Subcommittee, IEEE CIS Social Media Subcommittee
and IEEE SMC Technical Committee on Grey Systems.
---|---
|  Herman Ramon graduated as an agricultural engineer from Gent University. In
1993 he obtained a Ph.D. in applied biological sciences at the Katholieke
Universiteit Leuven. He is currently Professor at the Faculty of Agricultural
and Applied Biological Sciences of the Katholieke Universiteit Leuven,
lecturing on agricultural machinery and mechatronic systems for agricultural
machinery. He has a strong research interest in precision technologies and
advanced mechatronic systems for processes involved in the production chain of
food and nonfood materials, from the field to the end user.He is author or co-
author of more than 40 papers.
---|---
|  Wouter Saeys is currently Assistant Professor in Biosystems Engineering at
the Department of Biosystems of the University of Leuven in Belgium. He
obtained his Ph.D. at the same institute and was a visiting postdoc at the
School for Chemical Engineering and Advanced Materials of the University of
Newcastle upon Tyne, UK and at the Norwegian Food Research Institute - Nofima
Mat in Norway. His main research interests are optical sensing, process
monitoring and control with applications in food and agriculture. He is author
of 50 articles (ISI) and member of the editorial board of Biosystems
Engineering.
---|---


∎

11institutetext:  1The von Karman Institute for Fluid Dynamics, EA Department,
Waterloosesteenweg 72, Sint-Genesius-Rode 1640, Belgium
2Université libre de Bruxelles, Transfers, Interfaces and Processes (TIPs),
Av. Franklin Roosevelt 50, Brussels 1050, Belgium

# Multi-objective optimization of the magnetic wiping process in dip-coating

Fabio Pino1,2 Benoit Scheid2 Miguel A. Mendez1 Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

Electromagnetic wiping systems allow to pre-meter the coating thickness of the
liquid metal on a moving substrate. These systems have the potential to
provide more uniform coating and significantly higher production rates
compared to pneumatic wiping, but they require substantially larger amounts of
energy. This work presents a multi-objective optimization accounting for (1)
maximal wiping efficiency (2) maximal smoothness of the wiping meniscus, and
(3) minimal Joule heating. We present the Pareto front, identifying the best
wiping conditions given a set of weights for the three competing objectives.
The optimization was based on a 1D steady-state integral model, whose
prediction scales according to the Hartmann number (Ha). The optimization uses
a multi-gradient approach, with gradients computed with a combination of
finite differences and variational methods. The results show that the wiping
efficiency depends solely on Ha and not the magnetic field distribution.
Moreover, we show that the liquid thickness becomes insensitive to the
intensity of the magnetic field above a certain threshold and that the current
distribution (hence the Joule heating) is mildly affected by the magnetic
field’s intensity and shape.

##  1 Introduction

Hot-dip galvanization consists in coating a steel substrate with a thin layer
of zinc to prevent corrosion (Kistler and Schweizer, 1997). The substrate
usually takes the form of a large strip, dipped and then withdrawn from a bath
of molten zinc. The liquid excess is removed via a wiping system. In
galvannealing, the strip is heated for a few seconds after the wiping region
in an induction furnace (annealing furnace) when the zinc is still liquid.
Thereby, the iron in the steel and the zinc coating diffuse into each other,
forming a zinc-iron inter-metallic compound with a composition of
approximately $90\%$ iron and $10\%$ zinc in the bulk (Goodwin, 2013; Center,
2014).

The wiping is usually performed through impinging gas jets in a process called
jet wiping (Gosset and Buchlin, 2006; Hocking et al., 2011). The process is
simple and economically attractive, but it develops hydrodynamic instabilities
when the substrate speed exceeds 2 or 3 m/s, depending on the wiping strength.
These instabilities include splashing of the wiped liquid (Gosset and Buchlin,
2006) and oscillations of the impinging jet, which results in coating
undulations (Gosset et al., 2019; Mendez et al., 2019; Barreiro-Villaverde et
al., 2021).

An alternative to jet wiping is electromagnetic wiping, extensively
investigated by Bregand et al. (2010). These systems are based on the
ElectroMagnetic Brake (EMB) principle, that is the use of electromagnetic
forces to control the thickness of the liquid coating (Ernst et al., 2009;
Lloyd-Jones, C and Barker, HA and Worner, 1998; Dumont et al., 2011). This is
achieved using a strong magnetic field produced by a permanent magnet or high-
frequency induction generators. The relative motion of the liquid film and the
magnetic field create an induced current in the liquid film, which generates a
downward Lorentz force that wipes the excess of liquid. These systems have
been considered as “pre-wiping” devices in lines that target very large
substrate speeds (Lloyd-Jones, C and Barker, HA and Worner, 1998). Despite
offering promising performances and coating uniformity at large production
rates, the major challenges in their large-scale deployment are the energy
requirements and the possible overheating of the liquid zinc, which increases
the risks of promoting iron-zinc phase growth. The design of an
electromagnetic system must thus strike a balance between the contrasting
objectives of having sufficient “wiping capabilities” while minimizing
overheating and energy consumption.

This work investigates the trade-off between these competing goals using
multi-objective optimization. Specifically, we present the Pareto front,
defining the optimum according to three criteria: (1) wiping capabilities, (2)
Joule heating and (3) smoothness of the wiping liquid meniscus. The
optimization is based on a simple steady-state model, neglecting surface
tension. The remaining of the article is structured as follows. Section 4
presents the simplified model used in the optimization; this is based on the
steady-state streamwise momentum balance of the liquid film with the magnetic
effect expressed by the Hartmann number. Section 5 reports on the methods for
computing the relevant first-order sensitivity of the problem, which combined
variational methods and finite differences depending on the size of the
parameter vector. Section 6 defines the constrained multi-objective
optimization problem and the solution method based on a multi-gradient
descent. Finally, conclusions and perspectives are reported in Section 8.

Figure 1: Schematic of the hot-dip galvannealed coating process with the metal
strip dip into a bath of molten zinc, withdrawn vertically at a constant
velocity $U_{p}$. Outside the bath, the liquid excess is wiped out by a
magnetic actuator, and then the strip is heated in a galvannealing furnace.

##  2 Problem description

We consider liquid zinc with density $\rho=6570$kg/m3, dynamic viscosity
$\mu=3.5\times 10^{-3}$ Pa$\cdot$s, kinematic viscosity $\nu=\mu/\rho$,
surface tension $\sigma=0.7$N/m and electrical conductivity
$\sigma_{M}=0.028\times 10^{8}$ S/m. These properties are taken at around 500
Celsius degrees as in Giordanengo et al. (2000). Figure 1 shows the schematic
of the process: a metal strip is withdrawn from a bath of molten zinc at a
constant speed $U_{p}$ under the effect of gravity $g$, wiped by permanent
magnets (light blue circle) and then heated in a galvannealing furnace (green
box).

The reference system is centered on the substrate, with the positive
streamwise axis $x$ pointing downwards, the wall-normal axis $y$ towards the
free surface, and the cross-stream axis $z$ in the spanwise direction. We
assume a two-dimensional velocity vector field within the liquid, given by
$\mathbf{V}=[u(x,y),v(x,y),0]^{T}$. The flow is bounded by the substrate at
$y=0$, and by the free surface $h$, with the streamwise flow rate $q$ defined
as:

$q=\int_{0}^{h}u(y)\;dy.$ (1)

The magnetic field generated by the actuators is given by the vector:

$\mathbf{B}=[0,-B(x),0]^{T},$ (2)

where $B(x)$ is always positive since, for simplicity, we are only considering
the negative pole. The relative motion of the liquid zinc with respect to the
external magnetic field generates an induced current in the bulk. Neglecting
the electric potential difference, this is given by the vector:

$\mathbf{J}=\mathbf{V}\times\mathbf{B}=[0,0,j_{z}]^{T}=[0,0,-\sigma_{M}uB],$
(3)

The integraction of the induced current with the external magnetic field
results in a downward Lorentz force $\mathbf{F}_{L}$ defined as the cross
product of $\mathbf{B}$ and $\mathbf{J}$:

$\mathbf{F}_{L}=\mathbf{J}\times\mathbf{B}=[F_{L},0,0]^{T},$ (4)

with $F_{L}=\sigma_{M}uB^{2}$.

##  3 Scaling quantities

To recast the governing equations in nondimensional form, we scale the flow
velocity by the substrate velocity $U_{p}$ and lengths by the liquid film
thickness resulting from the steady state viscous-gravity balance (Mendez et
al., 2021):

$[u]=U_{p}\qquad\qquad[h]=\sqrt{\frac{\nu U_{p}}{g}}\,,$ (5)

where the $[\cdot]$ denotes reference quantities. Based on $[u]$ and $[h]$, we
obtain the scale for the flow rate as:

$[q]=[u][h]=\sqrt{\frac{\nu U_{p}^{3}}{g}}$ (6)

We scale the streamwise coordinate $x$ assuming long-wave conditions:

$[x]=[h]/\epsilon,$ (7)

with $\epsilon\ll 1.$

We define the magnetic field and the induced current reference quantities with
their maximum values:

$[B]=sup(B(x))\qquad\qquad[J]=\sigma_{M}U_{p}[B],$ (8)

We define the power density $P_{D}$ reference quantity as:

$[P_{D}]=U_{p}\rho g.$ (9)

The non-dimensional group governing the system is the Hartman number:

$\mbox{\text{Ha}}=\sqrt{\frac{\sigma_{M}[B]^{2}U_{p}}{\rho g}},$ (10)

which weighs the importance of electromagnetic to gravitational forces. In the
following, we use the hat notation $\hat{a}=a/[a]$ for the non-dimensional
quantities and the shorthand notation $\partial_{x}h=\partial h/\partial_{x}$
for partial derivatives.

##  4 Simplified model

In line with classical simplified engineering models for the jet wiping
(Thornton and Graff, 1976; Tuck, 1983), we assume a gently sloped interface
($\partial_{x}h\sim 0$) (Mendez et al., 2021), a unidimensional flow
($v(x,y)=0$), without surface tension (Dumont et al., 2011) and without
pressure gradient along the wall-normal and the streamwise direction. The
governing equation is given by the steady-state streamwise momentum balance,
which in non-dimensional form reads:

$\partial_{\hat{y}\hat{y}}\hat{u}+1-\mbox{\text{Ha}}^{2}\hat{B}^{2}\hat{u}=0,$
(11)

where the first term is the viscous stress, the second the gravitational
forces, and the last the Lorentz force. The non-dimensional magnetic field
$\hat{B}$ is modelled as a Gaussian along the streamwise direction, centred in
$\hat{x}=0$ with standard deviation $\gamma$:

$\hat{B}=\exp\Big{(}-\frac{\hat{x}^{2}}{2\gamma^{2}}\Big{)}.$ (12)

We impose a no-slip boundary condition at the wall, and a shear-free boundary
condition at the free surface (Kalliadasis et al., 2011, Chapter 2):

$\hat{u}(\hat{y}=0)=-1,\qquad\qquad\partial_{\hat{y}}\hat{u}(\hat{y}=\hat{h})=0.$
(13)

The solutions of (11) with boundary conditions (13) is given by:

$\begin{split}\hat{u}=\left(1+\frac{1}{\mbox{\text{Ha}}^{2}\hat{B}^{2}}\right)&\Big{(}\tanh{(\mbox{\text{Ha}}\hat{h}\hat{B})}\sinh{(\mbox{\text{Ha}}\hat{y}\hat{B})}-\cosh{(\mbox{\text{Ha}}\hat{y}\hat{B})}\Big{)}+\frac{1}{\mbox{\text{Ha}}^{2}\hat{B}^{2}},\end{split}$
(14)

as shown by Dumont et al. (2011). Integrating the streamwise velocity over the
liquid film thickness, we obtain a relation between $\hat{q}$ and $\hat{h}$:

$\begin{split}\hat{q}=g(\hat{h};\mbox{\text{Ha}}\hat{B})=&\frac{1}{\mbox{\text{Ha}}^{3}\hat{B}^{3}}\Big{\\{}\mbox{\text{Ha}}\,,\hat{h}\hat{B}-\tanh{(\mbox{\text{Ha}}\,\hat{h}\hat{B})}\left(\mbox{\text{Ha}}^{2}\hat{B}^{2}+1\right)\Big{\\}}.\end{split}$
(15)

Figure 2: Relation between the nondimensional flow rate $\hat{q}$ and liquid
film thickness $\hat{h}$ for different value of the nondimensional magnetic
field $\hat{B}$.

For $\mbox{\text{Ha}}\rightarrow 0$, (14) and (15) retrieve the relations
without wiping (Mendez et al., 2021). Figure 2 plots (15) for a fixed Ha at a
fixed $\hat{x}$ location for different values of the magnetic field $\hat{B}$.
In wiping conditions, $\hat{B}$ is varying along $\hat{x}$ with a minimum
$\hat{B}=0$ far from the actuator and a maximum $\hat{B}=1$ right at the pole
location ($\hat{x}=0$). The thickness $\hat{h}$ reaches a minimum at
$\hat{x}=0$. Its value is found by solving
$\partial_{\hat{h}}g(\hat{h};\mbox{\text{Ha}},\hat{B}=1)=0$, which yields:

$\hat{h}_{min}=\frac{1}{\mbox{\text{Ha}}}\operatorname{arccsc}\Big{(}\sqrt{\frac{1}{\mbox{\text{Ha}}^{2}+1}}\Big{)}.$
(16)

The flow rate $\hat{q}_{opt}$ is obtained injecting (16) in (15). The liquid
film profile $\hat{h}(\hat{x})$ is given by solving the equation:

$\hat{q}_{opt}-g(\hat{h}(\hat{x});\mbox{\text{Ha}},\hat{B}(\hat{x}))=0,$ (17)

at different locations along the stream-wise coordinate $\hat{x}$. In this
work, we consider $N=100$ equispaced points in a domain $\hat{x}\in[-10,10]$.

##  5 Sensitivity derivatives

The gradient-based multiobjective optimization carried out in this work
requires the computation of the jacobian (variation) of the state variables
$\hat{h}$, $\hat{q}$, $\hat{u}(\hat{y})$ and $\hat{j}$ with respect to the
magnetic wiping parameters $\mathbf{p}=[\mbox{\text{Ha}},\hat{B}]^{T}$:

$\partial_{\mathbf{p}}\hat{h}(\hat{x}),\qquad\qquad\partial_{\mathbf{p}}\hat{u}(\hat{y}),\qquad\qquad\partial_{\mathbf{p}}\hat{q},\qquad\qquad\partial_{\mathbf{p}}\hat{j}(\hat{x},\hat{y})\,.$
(18)

For a given Ha and $\hat{B}(x)$, the gradient $\partial_{\mathbf{p}}\hat{h}$
was computed with finite differences using the python library
scipy.optimize.approx_fprime on a set of $N$ equispaced points:

$\hat{h}(\hat{x}_{i})=\frac{\hat{h}(\hat{x}_{i};\mathbf{p}+\delta)-\hat{h}(\hat{x}_{i};\mathbf{p})}{\delta},$
(19)

with $\delta=10^{-8}$.

The total variation of $\hat{u}$ ($d\hat{u}/d\mathbf{p}$) is given by:

$\frac{d\hat{u}}{d\mathbf{p}}=\partial_{\mathbf{p}}\hat{u}+\partial_{\hat{h}}\hat{u}\partial_{\mathbf{p}}\hat{h}\,,$
(20)

where $\partial_{\hat{h}}\hat{u}$ was obtained differentiating (14):

$\partial_{\hat{h}}\hat{u}=\left(1+\frac{1}{\mbox{\text{Ha}}^{2}\hat{B}^{2}}\right)\frac{\mbox{\text{Ha}}\hat{B}}{\cosh{(\mbox{\text{Ha}}\hat{h}\hat{B})}^{2}}\sinh{(\mbox{\text{Ha}}\hat{y}\hat{B})}\,.$
(21)

The sensitivity gradient $\partial_{\mathbf{p}}\hat{u}(\hat{y})$ was
calculated using a variational method (Calver and Enright, 2017). This
consists in differentiating (11) with respect to $\mathbf{p}$, which leads to
a system for the sensitive gradients given by:

$\displaystyle\partial_{\mathbf{p}}\frac{d\mathbf{u}(\hat{y})}{d\hat{y}}\,=\,\mathbf{F}_{\mathbf{u}}\mathbf{S}_{\mathbf{p}}+\mathbf{F}_{\mathbf{p}}\,,$
(22)

where
$\mathbf{u}=[u_{1},u_{2}]^{T}=[\hat{u}(\hat{y}),\partial_{\hat{y}}\hat{u}(\hat{y})]^{T}$,
and $\mathbf{F}_{\mathbf{u}}$ and $\mathbf{F}_{\mathbf{p}}$ are the Jacobian
matrices with respect to $\mathbf{u}$ and $\mathbf{p}$ respectively:

$\mathbf{F_{u}}=\begin{pmatrix}0&1\\\
\mbox{\text{Ha}}^{2}\hat{B}^{2}&0\end{pmatrix}\,\qquad\mathbf{F}_{\mathbf{p}}=\begin{pmatrix}0&0\\\
2\mbox{\text{Ha}}\hat{B}^{2}u_{1}&2\mbox{\text{Ha}}^{2}\hat{B}u_{1}\end{pmatrix},$
(23)

and $\mathbf{S_{p}}$ is the matrix of sensitivity derivatives of $u_{1}$ and
$u_{2}$:

$\mathbf{S_{p}}=\begin{pmatrix}\partial_{\mbox{\text{Ha}}}u_{1}&\partial_{\hat{B}}u_{1}\\\
\partial_{\mbox{\text{Ha}}}u_{2}&\partial_{\hat{B}}u_{2}.\end{pmatrix}.$ (24)

Solving (11) and (22) with boundary conditions:

$\begin{split}u_{1}(\hat{y})|_{\hat{y}=0}=-1,\quad
u_{2}(\hat{y})|_{\hat{y}=\hat{h}},=0,\quad\partial_{\mbox{\text{Ha}}}u_{1}|_{\hat{y}=0}=0,\\\
\partial_{\hat{B}}u_{1}|_{\hat{y}=0}=0,\quad\partial_{\mbox{\text{Ha}}}u_{2}|_{\hat{y}=\hat{h}}=0,\quad\partial_{\hat{B}}u_{2}|_{\hat{y}=\hat{h}}=0.\end{split}$
(25)

gives the vector of sensitivity derivatives
$\partial_{\mathbf{p}}\hat{u}=[\partial_{\mbox{\text{Ha}}}\hat{u},\partial_{\hat{B}}\hat{u}]^{T}$.

The variation of the flow rate $\partial_{\mathbf{p}}\hat{q}$ is obtained as:

$\frac{d\mathbf{q}}{d\mathbf{p}}=\frac{d}{d\mathbf{p}}\int_{0}^{\mathbf{h}}\,\hat{u}(\hat{y}^{\prime})\,d\hat{y}^{\prime}=\hat{u}(\hat{h})\partial_{\mathbf{p}}\hat{h}+\int_{0}^{\hat{h}}\,\partial_{\mathbf{p}}\hat{u}(\hat{y}^{\prime})\,dy^{\prime}.$
(26)

Based on its definition in (3), the current density sensitivities are given
by:

$\partial_{\mbox{\text{Ha}}}\hat{j}=-\partial_{\mbox{\text{Ha}}}\hat{u}(\hat{y})\hat{B}(\hat{x}),\qquad\qquad\partial_{\hat{B}}\hat{j}=-\hat{u}(\hat{y})-\partial_{\hat{B}}\hat{u}(\hat{y})\hat{B}(\hat{x}).$
(27)

##  6 Multiobjectives optimization

We set up a constrained multiobjectives minimization problem problem for the
parameters Ha and $\gamma$, controlling the intensity and the shape of the
magnetic field distribution. The constraints impose that
$\mbox{\text{Ha}}<10$, above which the magnetic field becomes generally overly
large, a positive standard deviation $\gamma>0$, and a negative flow rate
$\hat{q}<0$ (that is some liquid is withdrawn from the plate):

$\boxed{\begin{aligned}
\min_{\mbox{\text{Ha}},\gamma}\quad&J(\mbox{\text{Ha}},\gamma)\\\
\textrm{s.t.}&\gamma>0,\,\,\mbox{\text{Ha}}<10\,\,\hat{q}<0\\\ \end{aligned}}$
(28)

The cost function $J$ adds up the contribution of three terms, weighting the
wiping strength ($J_{1}$), the smoothness of the free surface ($J_{2}$) and
the Joule heating due to the induced currents ($J_{3}$):

$J(Ha,\gamma)=J_{1}+J_{2}+J_{3}\,,$ (29)

where:

$\begin{split}J_{1}=-\hat{q},\quad
J_{2}=\frac{1}{2}\int\limits_{x_{1}}\limits^{x_{2}}\Big{(}\partial_{\hat{x}\hat{x}}\hat{h}(\hat{x})\Big{)}^{2}\;d\hat{x},\quad
J_{3}=\frac{1}{2}\int\limits_{x_{1}}\limits^{x_{2}}\int\limits_{0}\limits^{\hat{h}}\mbox{\text{Ha}}^{2}\hat{j}(\hat{x},\hat{y})^{2}\;d\hat{y}d\hat{x}\,.\end{split}$
(30)

Reducing the term $J_{1}$ requires reducing the flow rate and hence the liquid
film thickness. Reducing the term $J_{2}$ requires reducing the curvature of
the wiping meniscus, hence producing smoother transition from the thin film
solution at $x\rightarrow-\infty$ and the thick solution at
$x\rightarrow\infty$. A smoother transition implies smoother curvature of the
flow streamlines within the film, smaller velocity gradients and thus a more
stable film. The second derivative was computed by differentiating a fourth-
order spline, which interpolated the solution in all the grid points.

The term $J_{3}$ is related to the Joule heating due to the induced current
$\mathbf{J}$. In its differential form, the heating power $P_{D}$ is equal to
the dot product of the current density vector $\mathbf{{J}}$ times the
electric field $\mathbf{E}$. Since the current density and the electric field
are linked by the relation $\mathbf{E}=\mathbf{J}\sigma_{M}$, this yields the
following expression for $P_{D}$:

$P_{D}=\mathbf{J}\cdot\mathbf{E}=\mathbf{J}\cdot\mathbf{J}/\sigma_{M}=||\mathbf{J}||^{2}/\sigma_{M}\,,$
(31)

which in non-dimensional form reads (see also Smolentsev et al. (1997) and Mao
et al. (2008)):

$\hat{P}_{D}=\frac{\sigma_{M}U_{p}^{2}[B]^{2}}{[P_{D}]}\hat{j}^{2}=\mbox{\text{Ha}}^{2}\hat{j}^{2}.$
(32)

###  6.1 Multi-Gradient Descent method

The constrained optimization problem (28) is solved using the Stochastic
Multi-Subgradient Descent Algorithm (SMSGDA) method (Milojkovic et al., 2019).
This consists in updating the value of the parameters, starting from an
initial value $\mathbf{p}_{init}$, in a gradient descent fashion:

$\mathbf{p}^{\\{i+1\\}}=\mathbf{p}^{\\{i\\}}-\eta\nabla_{\mathbf{p}}J(\mathbf{p})$
(33)

where the superscript{i} indicates the number of optimization iterations,
$\eta$ is a small number controlling the parameter update and
$\nabla_{\mathbf{p}}J(\mathbf{p})$ is the common descent vector computed as
described in the following. This vector is given by the linear combination of
the gradient of the ith objective function
$\nabla_{\mathbf{p}}J_{i}(\mathbf{p})=[\partial_{\mbox{\text{Ha}}}J_{i},\partial_{\gamma}J_{i}]^{T}$
(see A for more details):

$\nabla_{\mathbf{p}}J(\mathbf{p})=\sum_{i=1}^{3}\,\alpha_{i}\nabla_{\mathbf{p}}J_{i}(\mathbf{p}),$
(34)

where $\alpha_{i}$ are the weighting coefficients. The values of these weights
are obtained by solving a quadratic optimization problem at each optimization
iteration:

$\min_{\alpha,\cdots,\alpha_{n}}\quad\Big{[}||\sum_{i=1}^{n}\,\alpha_{i}\nabla_{\mathbf{p}}J_{i}(\mathbf{p})||^{2}\;|\sum_{i=1}^{n}\,\alpha_{i}=1,\alpha_{i}\geq
0\Big{]}\\\ $ (35)

where $||\cdot||$ is the Euclidean norm. Due to the different scales of the
objective functions, the gradients can have significantly different norms
along the two objectives, hence causing these two to dominate on the other
one. To mitigate this risk, we normalize the gradients as follows:

$\nabla_{\mathbf{p}}\hat{J}_{i}(\mathbf{p})=\frac{\nabla_{\mathbf{p}}J_{i}(\mathbf{p})}{J_{i}(\mathbf{p}_{init})}\,,$
(36)

where $\nabla_{\mathbf{p}}\hat{J}_{i}(\mathbf{p})$ is the normalized gradient
vector, $\nabla_{\mathbf{p}}J_{i}(\mathbf{p})$ is the non-normalized gradient
and $J_{i}(\mathbf{p}_{init})$ is the initial loss related to the initial
parameter guess $\mathbf{p}_{init}$.

The optimization is concluded when the necessary optimality Karush-Kuhn-Tucker
(KKT) condition is met. This reads:

$\alpha,\cdots,\alpha_{n}\geq
0\qquad,\sum_{i=1}^{n}\,\alpha_{i}=1\qquad,\sum_{i=1}^{n}\,\alpha_{i}\nabla_{\mathbf{p}}J_{i}(\mathbf{p})=0\,.$

Every solution satisfying these conditions is called Pareto Stationary
(Désidéri, 2012). To find the Pareto optimal solutions, we run several
simulations starting from different $\mathbf{p}_{init}$ until they converged
on the point of the Pareto front.

##  7 Results

This section reports the sensitivity derivatives (Subsection 7.1) and the
Pareto front resulting from the multi-objective optimization problem in the
parameter $\mbox{\text{Ha}}-\gamma$ and the objective functions
$J_{1}-J_{2}-J_{3}$ spaces (Subsection 7.2).

###  7.1 Sensitivity derivatives

Figure 3: Examples of steady-state solution for different intensities of the
magnetic fields with the appearance of a hump, for different values of the
Hartman number Ha.

Figure 3 shows the liquid film for a magnetic Gaussian with $\gamma=3$ for
different Ha. As Ha increases, a portion of the film upstream of the actuator
curves, creating a local peak that becomes higher and steeper as Ha increases.
A hump starts to appear for Ha larger than two, approximately. The velocity
field beneath this hump is further analyzed in Subsection 7.2, where it is
illustrated for some of the optimal solutions belonging to the Pareto front.

####  7.1.1 Liquid film thickness

Figure 4 shows the sensitivity of the film thickness $\hat{h}$ with respect to
the Hartmann number Ha for two different magnetic field standard deviations:
$\gamma=1$ and $\gamma=3$. In both cases, we remark that $\hat{h}$ becomes
insensitive to the magnetic field for $\mbox{\text{Ha}}>4$, apart from a small
region of high sensitivity where a stationary hump starts to form. In
addition, we notice that as $\gamma$ grows, the hump start to form at a higher
value of Ha. The sensitivity presents a non-symmetric behaviour around the
symmetry axis of the magnetic field ($\hat{x}=0$). It has a larger negative
region that passes over the centreline of the magnetic Gaussian, affecting
points of the wiping region closer to the ran-back. A smaller $\gamma$ leads
to the same qualitative behaviour but with steeper gradients.

(a)

(b)

Figure 4: Colour map of the liquid film height sensitivity with respect to the
Hartmann number $Ha$ for a Gaussian magnetic field centred in zero (red dashed
line) with a standard deviation of (a) 1 and (b) 3. The horizontal axis shows
the one-dimensional streamwise coordinate $\hat{x}$, and the vertical axis
shows the different values of $H_{a}$ and the colours show the gradient.

####  7.1.2 Induced current density

Figure 5 shows the sensitivity of the electric current density function
$j(\hat{x},\hat{y})$ with respect to the Hartmann number at (a and b) $Ha=3$
and (c and d) $Ha=10$ and with respect to the magnetic Gaussian standard
deviations (a and c) $\gamma=3$ and (b and d) $\gamma=5$. For small Ha, the
most sensitive region is close to the peak of the magnetic Gaussian at the
middle of the liquid film. The region of high sensitivity extends around this
maximum, mainly in the run-back region of the film. It is interesting to note
that the most sensitive area in this region is close to the free surface. In
the far-field part, the sensitivity has a negative sign and is very small in
magnitude. As for the thickness sensitivity, an increase of Ha makes the
current density’s sensitivity very small in the whole domain apart from a
small region close to the magnetic peak and in the run-back.

(a)

(b)

(c)

(d)

Figure 5: Colour map of the induced current sensitivity sensitivities for
Hartmann number (a and b) $Ha=3$ and (c and d) $\mbox{\text{Ha}}=10$ and for
magnetic field standard deviations (a and c) $\gamma=3$ and (b and d)
$\gamma=5$.

###  7.2 Multi-objectives optimization

Figure 6 show $J_{1}$, $J_{2}$ and $J_{3}$ as a function of $Ha$ and $\gamma$.
We notice their evident conflicting behavior from their different shapes. The
curvature term $J_{2}$ tends to increase with Ha and decrease with $\gamma$,
peaking at $\mbox{\text{Ha}}=10$. The joule heating term $J_{3}$ grows with Ha
and $\gamma$. As seen in Subsec.4, the flow rate depends solely on Ha with a
monotone behavior that tends to become less steep for $\mbox{\text{Ha}}>4$.

(a)

(b)

(c)

Figure 6: Plot of the non-dimensional flow rate $\hat{q}$ the Joule heating
(b) and of the liquid film curvature (c) as a function of the Hartmann number
$H_{a}$ and the magnetic field standard deviation $\gamma$

By running the multi-objective optimization from different initial conditions,
we found the optimal parameters respecting the Karush-Kuhn-Tucker optimality
condition. Figure 7 shows the Pareto front collecting the optimal solutions in
the objective space (a) and the parameter space (b). The markers with the same
red shad in the two figures correspond to the same optimal solution. What
stands out from these plots is the role of $\gamma$. For $\mbox{\text{Ha}}<2$,
an increase of Ha is accompanied by an increase of $\gamma$. A larger magnetic
Gaussian compensates for the steepening of the interface due to Ha, while at
the same time keeping the Joule effect to a reasonably low value. For
$\mbox{\text{Ha}}>2$, an increase of Ha is accompanied by a decrease of
$\gamma$. In this case, it is better to accept a steep film interface with a
hump but with limited Joule heating. For larger, Ha the optimal $\gamma$
approach a plate at $\gamma\approx 8$. The Pareto front in Figure 7 is
tabulated in 2 for the reader convenience.

(a)

(b)

Figure 7: Pareto’s front in the (a) objective space and (b) the parameter
space.

Figure 8 show a collection of solutions belonging to the Pareto front, with
parameter values reported in Table 1. As Ha increases, the hump starts to
appear. Inside the hump, the current changes sing becoming positive. This
means that the streamwise velocity changes sign in this region. This suggests
the presence of a recirculation region under the hump, restricting the area
where the flow can pass to the far-field part of the film. This creates an
intense induced current around the magnetic peak with a sharp gradient. As a
consequence, in this area and the far field, the liquid is strongly affected
by an intense Joule heating.

Table 1: Cost functions and parameters’ values of four optimal solutions belonging to the Pareto front. Point | $J_{1}$ | $J_{2}$ | $J_{3}$ | Ha | $\gamma$
---|---|---|---|---|---
(a) | 0.62 | 2.15 | 2$\times 10^{-6}$ | 0.53 | 16.61
(b) | 0.45 | 11.56 | 4$\times 10^{-6}$ | 1.46 | 23.18
(c) | 0.2 | 33.42 | 0.004993 | 4.5 | 7.44
(d) | 0.12 | 61.17 | 0.027693 | 7.91 | 7.29

(a)

(b)

(c)

(d)

Figure 8: Induced current colour plot of optimal solutions belonging to the
Pareto front.

Figure 9 shows the solutions that minimize the distance to the so-called ideal
point, the origin of the normalized objective space reference system. The
liquid film is smooth without any hump and with an intense but limited Joule
heating, especially in the proximity of the moving strip. However, considering
an operating velocity of $1[m/s]$, the dimensional final coating is at least
one order of magnitude larger than the one required.

Figure 9: Solution minimizing the sum of the three cost functions in terms of
film thickness and current distribution for the multi-objective optimization
problem.

##  8 Conclusion

In this work, we presented a multi-objective optimization for electromagnetic
wiping in galvannealing coating processing. The optimization balances the need
to reduce the film thickness while preserving a smooth wiping meniscus and
minimizing the Joule heating. The problem was solved via a multi-gradient
descent optimization, providing the Pareto front in the objectives’ and
parameters’ spaces.

What stands out from this analysis is that the liquid film overheating and the
steepness of the wiping meniscus increase for $Ha>3$. Depending on the coating
requirements, one can choose solutions favouring one objective over the other.
This result could pave the way to a hybrid solution where the jet wiping is
supported by magnetic wiping in the pre-wiping or a pre-heating phase before
the galvannealing furnace. This would represent a good compromise between
wiping efficiencies and an increase in the temperature of the liquid, which
would reduce the effort of the galvannealing furnace in creating the Fe-Zn
alloy.

##  Appendix A Appendix A: Computation of the cost functions gradients

To solve the optimization problem, we calculate the gradients of $J$ with
respect to the magnetic parameters. The partial derivative of $J$ with respect
to Ha is given by:

$\partial_{\mbox{\text{Ha}}}J=\partial_{\mbox{\text{Ha}}}J_{1}+\partial_{\mbox{\text{Ha}}}J_{2}+\partial_{\mbox{\text{Ha}}}J_{3}$
(37)

where $\partial_{\mbox{\text{Ha}}}J_{1}$ is easily obtained via (26):

$\partial_{\mbox{\text{Ha}}}J_{1}=\partial_{\mbox{\text{Ha}}}\hat{q}=\hat{u}(\hat{h})\partial_{\mbox{\text{Ha}}}\hat{h}+\int_{0}^{\hat{h}}\,\partial_{\mbox{\text{Ha}}}\hat{u}(\hat{y})\,d\hat{y}\,,$
(38)

The partial derivative $\partial_{\mbox{\text{Ha}}}J_{2}$ is given by the
expression:

$\partial_{\mbox{\text{Ha}}}J_{2}=\int_{x_{1}}^{x_{2}}\,\partial_{\hat{x}\hat{x}}\hat{h}\partial_{\hat{x}\hat{x}}(\partial_{\mbox{\text{Ha}}}\hat{h})\,d\hat{x},$
(39)

having considered a domain $x\in[x_{1},x_{2}]$

The partial derivative of the Joule heating term
$\partial_{\mbox{\text{Ha}}}J_{3}$ is computed using the current density
sensitivity $\partial_{\mbox{\text{Ha}}}\hat{j}_{z}$:

$\begin{split}\partial_{\mbox{\text{Ha}}}J_{3}=&\frac{1}{2}\int_{x_{1}}^{x_{2}}\,\mbox{\text{Ha}}^{2}\hat{j}(\hat{h})^{2}\partial_{\mbox{\text{Ha}}}\hat{h}\,d\hat{x}+\int_{x_{1}}^{x_{2}}\int_{0}^{\hat{h}}\,\mbox{\text{Ha}}\hat{j}(\hat{y})^{2}\,d\hat{y}d\hat{x}+\int_{x_{1}}^{x_{2}}\int_{0}^{\hat{h}}\mbox{\text{Ha}}^{2}\hat{j}\partial_{\mbox{\text{Ha}}}\hat{j}\,d\hat{x}.\end{split}$
(40)

The gradient of $J$ with respect to $\gamma$ is given by:

$\partial_{\gamma}J=\partial_{\gamma}J_{1}+\partial_{\gamma}J_{2}+\partial_{\gamma}J_{3},$
(41)

where:

$\partial_{\gamma}J_{1}=\partial_{\gamma}\hat{q}=\hat{u}(\hat{h})\partial_{\hat{B}}\hat{h}\partial_{\gamma}\hat{B}+\int_{0}^{\hat{h}}\,\partial_{\hat{B}}\hat{u}(\hat{y})\partial_{\gamma}\hat{B}\,d\hat{y}\,,$
(42a)
$\partial_{\gamma}J_{2}=\int_{x_{1}}^{x_{2}}\,\partial_{\hat{x}\hat{x}}\hat{h}\partial_{\hat{x}\hat{x}}(\partial_{\hat{B}}\hat{h})\partial_{\gamma}\hat{B}\,d\hat{x},$
(42b)
$\begin{split}\partial_{\gamma}J_{3}=&\frac{1}{2}\int_{x_{1}}^{x_{2}}\,\mbox{\text{Ha}}^{2}\hat{j}(\hat{h})^{2}\partial_{\hat{B}}\hat{h}\partial_{\gamma}\hat{B}\,d\hat{x}+\int_{x_{1}}^{x_{2}}\int_{0}^{\hat{h}}\,\mbox{\text{Ha}}^{2}\hat{j}\partial_{\hat{B}}\hat{j}\partial_{\gamma}\hat{B}\,d\hat{y}d\hat{x}.\end{split}$
(42c)

##  Appendix B Appendix B: Tabulated Values of the identified Pareto front

Table 2 tabulates the identified Pareto front shown in Figure 7.

Table 2: Value of the magnetic parameters ($Ha,\gamma$) and the objective functions ($J_{1},J_{2},J_{3}$) along the identified Pareto front. Ha | $\gamma$ | $J_{1}$ | $J_{2}$ | $J_{3}$
---|---|---|---|---
0.53, | 16.61, | 0.62, | 2e-06, | 2.15,
0.63, | 17.58, | 0.6, | 2e-06, | 2.9,
0.69, | 18.23, | 0.59, | 2e-06, | 3.43,
0.74, | 18.73, | 0.58, | 2e-06, | 3.96,
1.06, | 22.18, | 0.52, | 2e-06, | 7.21,
1.19, | 22.59, | 0.5, | 2e-06, | 8.59,
1.34, | 22.97, | 0.47, | 3e-06, | 10.23,
1.46, | 23.18, | 0.45, | 4e-06, | 11.56,
2.14, | 9.66, | 0.36, | 0.000261, | 15.61,
2.31, | 9.26, | 0.34, | 0.000368, | 16.88,
2.5, | 8.89, | 0.33, | 0.000518, | 18.23,
2.7, | 8.54, | 0.31, | 0.000733, | 19.68,
3.34, | 7.85, | 0.26, | 0.00174, | 24.34,
3.6, | 7.7, | 0.25, | 0.002283, | 26.25,
3.74, | 7.63, | 0.24, | 0.002646, | 27.38,
3.87, | 7.59, | 0.23, | 0.002965, | 28.38,
4.5, | 7.44, | 0.2, | 0.004993, | 33.42,
4.67, | 7.42, | 0.2, | 0.005651, | 34.81,
4.8, | 7.41, | 0.19, | 0.006165, | 35.87,
4.9, | 7.39, | 0.19, | 0.006634, | 36.7,
5.6, | 7.32, | 0.17, | 0.010218, | 42.26,
7.91 | 7.29 | 0.12 | 0.027693 | 61.17

###### Acknowledgements.

F.Pino is supported by an F.R.S.-FNRS FRIA grant, and his PhD research project
is funded by Arcelor-Mittal. B.Scheid is Research Director at F.R.S.-FNRS.

## References

  * Barreiro-Villaverde et al. (2021) Barreiro-Villaverde D, Gosset A, Mendez M (2021) On the dynamics of jet wiping: Numerical simulations and modal analysis. Physics of Fluids 33(6), DOI 10.1063/5.0051451
  * Bregand et al. (2010) Bregand O, Crahay J, Gerkens P (2010) New wiping techniques to efficiently produce suitable coating layers at high speed in the hot dip galvanising process. DOI 10.2777/86670
  * Calver and Enright (2017) Calver J, Enright W (2017) Numerical methods for computing sensitivities for odes and ddes. Numerical Algorithms 74(4):1101–1117
  * Center (2014) Center G (2014) 1.3 galvanneal – differences from galvanize. Tech. rep., GalvInfo Notes, International Zinc Association
  * Désidéri (2012) Désidéri JA (2012) Multiple-gradient descent algorithm (mgda) for multiobjective optimization. Comptes Rendus Mathematique 350(5-6):313–318
  * Dumont et al. (2011) Dumont M, Ernst R, Fautrelle Y, Grenier B, Hardy JJ, Anderhuber M (2011) New DC electromagnetic wiping system for hot-dip coating. COMPEL - The International Journal for Computation and Mathematics in Electrical and Electronic Engineering 30(5):1663–1671, DOI 10.1108/03321641111152810
  * Ernst et al. (2009) Ernst R, Fautrelle Y, Bianchi A, Iliescu M (2009) Working principle of an electromagnetic wiping system. Magnetohydrodynamics 45(1):59–71
  * Giordanengo et al. (2000) Giordanengo B, Moussa AB, Makradi A, Chaaba H, Gasser J (2000) New precise determination of the high temperature unusual temperature dependent thermopower of liquid divalent cadmium and zinc. Journal of Physics: Condensed Matter 12(15):3595
  * Goodwin (2013) Goodwin FE (2013) Developments in the production of galvannealed steel for automotive. Transactions of the Indian Institute of Metals 66(5):671–676
  * Gosset and Buchlin (2006) Gosset A, Buchlin JM (2006) Jet wiping in hot-dip galvanization. Journal of Fluids Engineering 129(4):466–475, DOI 10.1115/1.2436585
  * Gosset et al. (2019) Gosset A, Mendez MA, Buchlin JM (2019) An experimental analysis of the stability of the jet wiping process: Part I – Characterization of the coating uniformity. Experimental Thermal and Fluid Science 103:51–65, DOI 10.1016/j.expthermflusci.2018.12.029
  * Hocking et al. (2011) Hocking GC, Sweatman W, Fitt A, Breward C (2011) Deformations during jet-stripping in the galvanizing process. Journal of Engineering Mathematics 70:297–306
  * Kalliadasis et al. (2011) Kalliadasis S, Ruyer-Quil C, Scheid B, Velarde MG (2011) Falling liquid films, vol 176. Springer Science & Business Media
  * Kistler and Schweizer (1997) Kistler SF, Schweizer PM (1997) Coating science and technology: An overview. Liquid Film Coating pp 3–15
  * Lloyd-Jones, C and Barker, HA and Worner (1998) Lloyd-Jones, C and Barker, HA and Worner V (1998) Investigation into magnetic wiping techniques as alternative to gas wiping on hot dip galvanising lines. Ironmaking & steelmaking 25(2):117
  * Mao et al. (2008) Mao J, Aleksandrova S, Molokov S (2008) Joule heating in magnetohydrodynamic flows in channels with thin conducting walls. International journal of heat and mass transfer 51(17-18):4392–4399
  * Mendez et al. (2019) Mendez M, Gosset A, Buchlin JM (2019) Experimental analysis of the stability of the jet wiping process, part ii: Multiscale modal analysis of the gas jet-liquid film interaction. Experimental Thermal and Fluid Science 106:48–67
  * Mendez et al. (2021) Mendez M, Gosset A, Scheid B, Balabane M, Buchlin JM (2021) Dynamics of the jet wiping process via integral models. Journal of Fluid Mechanics 911
  * Milojkovic et al. (2019) Milojkovic N, Antognini D, Bergamin G, Faltings B, Musat C (2019) Multi-gradient descent for multi-objective recommender systems. arXiv preprint arXiv:200100846
  * Smolentsev et al. (1997) Smolentsev SY, Tananaev A, Dajeh D (1997) Laminar heat transfer in the mhd flow in a rectangular duct 1. numerical and analytical investigations. Magnetohydrodynamics c/c of Magnitnaya Gidrodinamika 33:148–154
  * Thornton and Graff (1976) Thornton JA, Graff HF (1976) An analytical description of the jet finishing process for hot-dip metallic coatings on strip. Metallurgical Transactions B 7(4):607–618
  * Tuck (1983) Tuck E (1983) Continuous coating with gravity and jet stripping. The Physics of fluids 26(9):2352–2358