# Virtual Reference Feedback Tuning with robustness constraints: a swarm
intelligence solution

Luan Vinícius Fiorio<EMAIL_ADDRESS>Chrystian Lenon Remes  Patrick Wheeler
Yales Rômulo de Novaes

###### Abstract

The simplified modeling of a complex system allied with a low-order controller
structure can lead to poor closed-loop performance and robustness. A feasible
solution is to avoid the necessity of a model by using data for the controller
design. The Virtual Reference Feedback Tuning (VRFT) is a data-driven design
method that only requires a single batch of data and solves a reference
tracking problem, although with no guarantee of robustness. In this work, the
inclusion of an $\mathcal{H}_{\infty}$ robustness constraint to the VRFT cost
function is addressed. The estimation of the $\mathcal{H}_{\infty}$ norm of
the sensitivity transfer function is extended to maintain the one-shot
characteristic of the VRFT. Swarm intelligence algorithms are used to solve
the non-convex cost function. The proposed method is applied in two real-world
inspired problems with four different swarm intelligence algorithms, which are
compared with each other through a Monte Carlo experiment of 50 executions.
The obtained results are satisfactory, achieving the desired robustness
criteria.

###### keywords:

Data-driven control , Robust control , Swarm intelligence algorithms , Virtual
Reference Feedback Tuning

[inst0]department=Department of Electrical Engineering,,
organization=Eindhoven University of Technology, postcode=5600 MB,
city=Eindhoven, country=the Netherlands

[inst1]department=Department of Electrical Engineering,, organization=Santa
Catarina State University, postcode=89219-710, city=Joinville, country=Brazil

[inst2]department=Department of Electrical Engineering,, organization=Federal
University of Paraná, postcode=80060-000, city=Curitiba, country=Brazil

[inst3]department=Department of Electrical Engineering,, organization=Federal
University of Rio Grande do Sul, postcode=90035-190, city=Porto Alegre,
country=Brazil

[inst4]department=Department of Electrical and Electronic Engineering,,
organization=University of Nottingham, postcode=NG7 2RD, city=Nottingham,
country=United Kingdom

##  1 Introduction

The complexity of certain processes usually requires simplification in the
mathematical modeling to allow proper controller design [1], e.g., for power
systems [2, 3]. For dc-dc converters, for example, the majority of the control
techniques assume the existence of an accurate model [4, 5], which can be
difficult to obtain since power converters have nonlinear dynamics. Low-order
controllers are commonly used in the industry, such as Proportional Integral
(PI) and Proportional Integral Derivative (PID) controllers [6, 7, 8, 9], but
a robust design can still be a challenge. Oversimplified modeling due to the
process complexity, and/or a limited performance of the chosen controller
structure are possible reasons for achieving poor performance and robustness
of the controlled system [10].

On the other hand, data-driven control design techniques are used to overcome
common problems related to system models, such as the dilemma on
representativity and complexity, or even the unavailability of those [11, 12,
13]. Some of the data-driven approaches require several plant experiments and
iterative acquisition of data, like Iterative Feedback Tuning (IFT) [14],
Correlation-based Tuning (CbT) [15], and more recent methods such as the deep
neural network-based implementation of data-driven Iterative Learning Control
(ILC) [16], and data-driven optimization-based ILC [17]. While other methods
like Virtual Reference Feedback Tuning (VRFT) [18], Direct Iterative Tuning
(DIT) [19], Optimal Controller Identification (OCI) [20], Virtual Disturbance
Feedback Tuning (VDFT) [21], Data-Driven Linear Quadratic Regulator (DD-LQR)
[22], and data-driven predictive model control [23] only require a single
batch of data (one-shot) in order to tune the controller parameters. If data
acquisition is difficult, its repetition is prevented by a tight deadline, or
memory constraints are an issue, one-shot solutions are usually preferred over
iterative solutions. For more details, a recent and systematic review of data-
driven control methods can be found at [24].

Robustness considering low-order controllers is a frequent topic of discussion
[25, 26] since specific processes can present uncertainties, as well as
disturbances that might occur over time. Another point to be observed is that
a poor choice of reference model or limited controller class in, e.g., the
VRFT design, may result in poor performance or robustness [27]. Since
robustness can be measured by the $\mathcal{H}_{\infty}$ norm of the
sensitivity transfer function $S(z)$ of a closed-loop system [28], its
inclusion in the data-driven design of controllers can be considered, allowing
for a more robust design when necessary. A recent methodology proposed in [29]
has suggested the inclusion of a robustness criterion in the VRFT design, at
the expense of: i) more experiments since the proposed design procedure
iterates in a trial-and-error fashion until the desired robustness is achieved
and, essentially removes one of the greatest advantages of the VRFT - being a
one-shot method; and ii) this type of iterative procedure requiring more
background knowledge from the designer for choosing reference models and
specifying requirements.

Also, a data-driven one-shot approach for multivariable systems regarding
robust solutions of $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$ criteria, and
loop-shaping specifications, has been recently addressed [30]. But the method
relies on an initial solution that influences the final one since it is
convexified by linearization around the initial stabilizing controller, which
is undesirable. In the case of data-driven one-shot approaches for Single-
Input Single-Output (SISO) systems, considering an $\mathcal{H}_{\infty}$
robust performance criteria, a recent work [31] proposes a technique where the
order of the controller is increased as the solution converges and only deals
with noise-free data. Both methods [30, 31] use only frequency-domain data.
Other data-driven robust solutions are achieved in an online fashion and
require higher computational processing than offline techniques, and also
present the demand of measuring signals in real-time. Some of those techniques
are the use of the modified Riccati equation with online data-driven learning
[32] and the application of a data-driven Model Predictive Control method with
robustness guarantees [33].

Among the aforementioned data-driven design techniques, the VRFT has been more
broadly applied to several classes of problems, a fact that can be useful to
attest to the feasibility of the proposed method of this paper, as it requires
less computational cost than, e.g., OCI and online methods. Therefore, this
work is based on the VRFT technique, as presented in [18] and [34], proposing
the inclusion of a robustness constraint, in terms of the
$\mathcal{H}_{\infty}$ norm of the sensitivity transfer function of the
system, in the VRFT design, but maintaining one of its most attractive
features: the necessity of only one batch of data. The constraint is included
in the VRFT cost function as a penalty [35], which compromises its convex
behavior. To deal with the non-convex optimization procedure, the proposed
method is tackled in two main steps: i) the design, in a data-driven fashion,
of a controller using the VRFT approach, if a previous controller is
unexistent; and ii) the application of a metaheuristic optimization algorithm
to minimize the cost function considering an $||S(z)||_{\infty}$ norm
constraint, using the controller from the previous step as an initial
solution.

Notice that metaheuristics, more specifically the Particle Swarm Optimization
(PSO) algorithm, have been used on the VRFT design for enhancing the
performance of the controller in terms of reference tracking on hybrid systems
[36]. Also, more closely related to this work, PSO has been used for solving
the VRFT with a modified cost function which includes the optimization of the
reference model based on a reference tracking metric [37]. The work focuses on
reducing the mismatch between the reference model and the actual controlled
system but does not mention robustness.

Since metaheuristic algorithms may work well for a certain class of problems,
but might fail over other problems, according to the No Free Lunch (NFL)
theorems [38], more than a single metaheuristic optimization algorithm should
be evaluated. Looking over the available types of metaheuristics, three can be
highlighted: evolutionary algorithms [39]; physics-based algorithms [40]; and
swarm intelligence algorithms. Although some authors group evolutionary
algorithms with multiple agents and swarm algorithms together [41], this paper
considers the two classes separately, as done by other authors in the
metaheuristic optimization subject [42]. This work focuses on swarm
intelligence algorithms since they usually have fewer parameters to be tuned
by the user or designer [41]. Four swarm intelligence algorithms are
considered: Particle Swarm Optimization [43] and Artificial Bee Colony (ABC)
[44] since those are two of the most used in literature; and the two swarm
intelligence algorithms with the least number of hyperparameters, Grey Wolf
Optimizer (GWO) [42], and its most recent version, the Improved Grey Wolf
Optimizer (I-GWO) [45].

In summary, the main contributions of this work are:

  * 1.

in theoretical terms, the derivation of the expressions of the signals for the
data-driven estimation of the impulse response of the sensitivity transfer
function of a system, in one-shot, at Subsection 4.1.1;

  * 2.

the proposal of a modified cost function to the VRFT method regarding the
inclusion of a robustness constraint, at Section 5, and the use of swarm-
intelligence algorithms for solving the proposed cost function;

  * 3.

at Section 6, the use and comparison of four different swarm-intelligence
algorithms - PSO, ABC, GWO, and I-GWO - in the proposed method for solving two
real-world inspired problems based on the structure of widely used DC-DC
converters.

This paper is structured as follows: Section 2 describes the system that is
considered in this paper for the theoretical formulation; Section 3 presents
the basis of data-driven controller design and details the basic VRFT design
procedure; Section 4 explains the method for the $\mathcal{H}_{\infty}$ norm
estimation of the sensitivity transfer function; Section 5 presents the
proposed method and details the four considered swarm intelligence algorithms;
Section 6 illustrates and validates the method by showing its application in
two real-world inspired examples; and finally, Section 7 concludes this work.

##  2 Preliminaries: description of the system

The system considered in this paper for the theoretical formulation is a
discrete-time, causal, linear time-invariant, and Single-Input Single-Output
(SISO) system $G(z)$. It is considered that $z$ is the forward discrete time-
shift operator such that $zx(k)=x(k+1)$. The output $y(k)$ of the system can
be described as

$y(k)=G(z)u(k)+v(k),$ (1)

where $u(k)$ is the input signal and $v(k)$ is the process noise - stochastic
effects that are not represented by $G(z)$, i.e., not captured by the input-
output relation of $u(k)$ and $y(k)$.

The closed-loop system taken into account in this work is regarded as a
controller $C(z)$ with the process $G(z)$ and a unit gain feedback, as shown
in Figure 1, where $r(k)$ is the reference signal and $e(k)$ is the error
signal. The closed-loop control law is

$u(k)=C(z,\rho)(r(k)-y(k)),$ (2)

with

$C(z,\rho)=\rho^{\prime}\bar{C}(z)$ (3)

being a controller with parameter $\rho\in\mathbb{R}^{p}$, $\bar{C}(z)$ a
vector of transfer functions belonging to the controller class $\mathcal{C}$
(e.g., PI or PID controller classes) and $r(k)$ the reference signal.

The output of the closed-loop system is given as

$y(k)=T(z)r(k)+S(z)v(k),$ (4)

where the reference signal $r(k)$ is applied to the transfer function from the
reference $r(k)$ to the output $y(k)$, $T(z)$, with

$T(z)=\frac{C(z)G(z)}{1+C(z)G(z)},$ (5)

and $S(z)$ is the sensitivity transfer function such that $S(z)+T(z)=1$ and

$S(z)=\frac{1}{1+C(z)G(z)}.$ (6)

$C(z)$$G(z)$$u(k)$$r(k)\ \ +$$v(k)$$+$$e(k)$$+$$y(k)$$-$ Figure 1: Block
diagram of the considered closed-loop system structure for this paper.

In the next section, the Model Reference Control (MRC) is introduced. The VRFT
method is described for the cases where the process has minimum or non-minimum
phase.

##  3 Virtual Reference Feedback Tuning (VRFT)

The Model Reference Control (MRC), which is the basis for the VRFT, more
generally called model matching control [46], concerns the problem of
reference tracking of the closed-loop system’s response, disregarding the
effects of noise at the output [27].

In order to obtain a controller, the MRC requires the designer to elaborate a
target transfer function for the controlled closed-loop system, called a
reference model ($T_{d}(z)$), which generates the output
$y_{d}(k)=T_{d}(z)r(k)$. A reference tracking performance criterion evaluated
by the two-norm tracking error is then obtained by solving the optimization
problem —l— ρJ^MR(ρ) = —— (T(z,ρ) - T_d(z)) r(k) ——_2^2, which can be solved
by considering (5), resulting in the solution controller for the MRC, called
the ideal controller $C_{d}(z)$.

The VRFT is a one-shot optimization data-driven controller design technique
based on the MRC. It is defined as one-shot since only a single batch of
input-output data is required to solve the model reference control problem
(3), which can be done by the use of least squares when the controller is
linearly parametrized as in (3), resulting in the parameter $\rho$ of a
controller with predefined class. The VRFT design depicted in this paper
follows the procedures of [27, 11].

Consider an experiment, in open-loop or closed-loop, that results in a batch
of collected data $\\{u,y\\}_{k=1}^{N}$. A virtual reference signal
$\bar{r}(k)$ is defined such that $T_{d}(z)\bar{r}(k)=y(k)$. A virtual error
can be obtained as $\bar{e}(k)=\bar{r}(k)-y(k)=(T_{d}^{-1}(z)-1)y(k)$. In
summary, a controller $C(z,\rho)=\rho^{\prime}\bar{C}(z)$ is considered
satisfactory if it generates $u(k)$ when fed by $\bar{e}(k)$. The closed-loop
of such block diagram for the VRFT controller design is illustrated in Figure
2.

Figure 2: Closed-loop block diagram for the VRFT controller design.

The VRFT solves the optimization problem —l— ρJ^VR(ρ) = —— u(k) - C(z,ρ)
(T_d^-1(z) - 1) y(k) ——_2^2, having the same minimum as (3) if the ideal
controller $C_{d}(z)$ in (5) belongs to the same controller class
$\mathcal{C}=\\{C(z,\rho),\rho\in\mathbb{R}^{p}\\}$ as $C(z,\rho)$. To
compensate for the fact that the ideal controller rarely belongs to the chosen
controller class, a filter $L(z)$ is applied to the data to approximate the
minimum of $J^{VR}$ to the minimum of $J^{MR}$, where the amplitude should
satisfy [27]:

$|L(e^{j\Omega})|^{2}=|T_{d}(e^{j\Omega})|^{2}|1-T_{d}(e^{j\Omega})|^{2}\frac{\Phi_{r}(e^{j\Omega})}{\Phi_{u}(e^{j\Omega})},\quad\forall\Omega\in[-\pi,\pi],$
(7)

where $x(e^{j\Omega})$, with $x$ representing any signal or system, represents
the Discrete Fourier Transform of $x(k)$,
$\Phi_{r}(e^{j\Omega}),\Phi_{u}(e^{j\Omega})$ are, respectively, the power
spectra of the signals $r(k),u(k)$.

Instrumental variables can be used in order to suppress estimation bias caused
by the noise in data [47], requiring the use of a second data batch. In
practice, the input signal can be formed by two identical sequences in the
same experiment, if memory restrictions do not impose a problem. Then, the
signals can be synced together afterward, resulting in two batches of data
from one single experiment.

In the presence of a Non-Minimum Phase (NMP) zero at the process, a flexible
reference model can be used, as presented in [27]. The optimization problem
(2), then, becomes —l— ρJ^VR(ρ) = —— η’ F(z) (u(k) + ρ’ ¯C(z) y(k)) - ρ’ ¯C(z)
y(k) ——_2^2, where $\eta\in\mathbb{R}^{m}$, and $F(z)$ is a vector of transfer
functions such that $T_{d}(z,\eta)=\eta^{\prime}F(z)$. The step-by-step design
for the VRFT with a flexible reference model, from data collection to the
algorithm design, is detailed in [11].

To be able to include a robustness criterion at the VRFT cost function, a
means of evaluating this parameter is required. Nonetheless, this is
approached in the next section.

##  4 Robustness index estimation

Depending on the choice of $T_{d}(z)$ or the controller class $\mathcal{C}$,
as well as the response of the plant $G(z)$, the VRFT-designed controller can
result in poor robustness of the controlled process. For such cases, a
robustness constraint can be included in the form of the
$\mathcal{H}_{\infty}$ norm of $S$, also called maximum sensitivity ($M_{S}$),
which can be used as a measure of robustness [28].

Typically, a system that presents $M_{S}>2$ is considered to have poor
robustness [28]. In this context, and considering the use of data-driven
design approaches, it is necessary to estimate the value of
$||S(z)||_{\infty}$ in a data-driven way since it is assumed that no plant
model is available to the designer. By this reasoning, the estimation of
$M_{S}$ is explained in the following subsection.

###  4.1 Estimation of $M_{S}$

The $\mathcal{H}_{\infty}$ norm estimation procedure developed in this work is
based on the Impulse Response (IR) of the system, as presented in [48],
modified from [49] to a SISO impulse response identification procedure, which
allows for a regularized estimation according to the existing literature [50].
Also, in order to maintain the one-shot characteristic of the VRFT, the
estimation of the $\mathcal{H}_{\infty}$ norm of $S$ based on impulse response
is addressed as follows.

Consider the linear discrete-time causal and SISO system $S$, represented by
its transfer function $S(z,\rho)$, such that its output signal $\psi(k)$ is
given by the equation

$S:\psi(k)=s(k)\ast\zeta(k)=\sum_{n=0}^{\infty}s(k-n)\zeta(n),$ (8)

where $\zeta(k)$ is the input signal of $S$, whose impulse response is $s(k)$.

Since (8) requires infinite data to be obtained, an order $M$ is defined such
that it is assumed that any IR term greater than $M$ is negligible, which is
valid for stable systems since $\lim_{k\to\infty}s(k)=0$. Nevertheless,
equation (8) can be truncated to $M$ terms, leading to:

$S:\psi(k)=\sum_{n=0}^{\infty}s(k-n)\zeta(n)\approx\underbrace{\sum_{n=0}^{M}s(k-n)\zeta(n)}_{|s(M+1)|<\epsilon\text{,
with }\epsilon\to 0^{+}}.$ (9)

On the other hand, the definition of the $\mathcal{H}_{\infty}$ norm, when
applied to the system $S$, can be written as [28]

$\mathcal{H}_{\infty}:||S||_{\infty}=\max_{\zeta(k)\neq
0}\dfrac{||s(k)\ast\zeta(k)||_{2}}{||\zeta(k)||_{2}},$ (10)

which requires the whole set of possible inputs $\\{\zeta(k)\neq 0\\}$.
Therefore, expression (10) cannot be directly calculated. An alternative
strategy is to obtain a matrix relation for $S$, which allows for the use of
induced norm properties.

Expanding (9) to the $M$ first terms gives

$\begin{cases}&\psi(0)=s(0)\zeta(0)\\\ &\psi(1)=s(1)\zeta(0)+s(0)\zeta(1)\\\
&\vdots\\\ &\psi(M)=s(M)\zeta(0)+\cdots+s(0)\zeta(M),\end{cases}$ (11)

and the following matrix relation truncated at $M$ elements is obtained:

$\underbrace{\begin{bmatrix}\psi(0)\\\ \psi(1)\\\ \cdots\\\
\psi(M)\end{bmatrix}}_{\Psi_{M}}=\underbrace{\begin{bmatrix}s(0)&0&\cdots&0\\\
s(1)&s(0)&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
s(M)&s(M-1)&\cdots&s(0)\end{bmatrix}}_{S_{M}}\underbrace{\begin{bmatrix}\zeta(0)\\\
\zeta(1)\\\ \cdots\\\ \zeta(M)\end{bmatrix}}_{Z_{M}}.$ (12)

From the assumption that the order $M$ is sufficiently high, it can be said
that matrix $S_{M}$ characterizes the IR $s(k)$ of $S$.

A useful matrix property is the induced norm [28, A.5], which can be applied
to (12), such that

$||S_{M}||_{ip}=\max_{Z_{M}\neq 0}\dfrac{||S_{M}Z_{M}||_{p}}{||Z_{M}||_{p}},$
(13)

where the subscript $i$ stands for induced. In short, (13) is a matrix form of
representing the system gain considering a set of possible input signals
$Z_{M}$. From the induced-2 norm, the following is obtained

$||S_{M}||_{i2}=\bar{\sigma}(S_{M})=\sqrt{\lambda_{max}\left(S_{M}^{\prime}S_{M}\right)},$
(14)

where $\bar{\sigma}$ and $\lambda_{max}$ stands for largest singular value and
largest eigenvalue, respectively, and comparing (10) with (13), it can be seen
that

$||S||_{\infty}\approx\max_{Z_{M}\neq
0}\frac{||S_{M}Z_{M}||_{2}}{||Z_{M}||_{2}}=\sqrt{\lambda_{\max}(S_{M}^{\prime}S_{M})}.$
(15)

Since the $||S(z,\rho)||_{\infty}$ norm can be estimated based on its IR via
(15), an expression for the input signal $\zeta(k)$ and the output signal
$\psi(k)$ of $S(z,\rho)$ has to be derived in order to estimate its impulse
response in the VRFT design context.

####  4.1.1 Input-output signals of the sensitivity transfer function

Considering the system presented in Figure 2, its sensitivity transfer
function in (6) can be rewritten as

$1+C(z,\rho)G(z)=S^{-1}(z,\rho).$ (16)

Assuming that $u(k)$ is sufficiently informative to capture all relevant
characteristics of $S(z,\rho)$, multiplying both sides of (16) by $u(k)$
achieves

$u(k)+C(z,\rho)G(z)u(k)=S^{-1}(z,\rho)u(k).$ (17)

It is known that $G(z)u(k)=y(k)$. Substituting such relation in (17):

$u(k)+C(z,\rho)y(k)=S^{-1}(z,\rho)u(k).$ (18)

From (18), the signals

$\xi(k)=u(k),\quad\zeta(k)=u(k)+C(z,\rho)y(k),$ (19)

can be derived. Expressions (18) and (19) mean that when a signal $\zeta(k)$
formed by $u(k)+C(z,\rho)y(k)$ is applied to $S(z,\rho)$, an output
$\xi(k)=u(k)$ is obtained. Therefore, the impulse response of $S(z,\rho)$ can
be estimated considering the data set $\\{\xi,\zeta\\}_{k=1}^{N}$.

In this work, the IR estimation is made through identification with
regularization techniques since: i) the variance of the estimates increases
with the order $M$, which is suppressed with the use of regularization [50];
and ii) knowing that IR is a sparse signal for sufficiently high $M$, the use
of regularization is known to improve sparse signal estimates [51]. The
algorithm for regularized estimation of the impulse response is described in
[50, 52], and is available in MATLAB® [53], Python [54], and R [55].

The inclusion of the $\mathcal{H}_{\infty}$ constraint in the cost function
spoils the convexity characteristic of the VRFT and the solution can no longer
be obtained through the least squares algorithm. A strategy to deal with local
minima and other characteristics that may arise from a non-convex cost
function is to use metaheuristic optimization [56]. Nevertheless, the VRFT
with robustness constraints is formally proposed in the next section, which
also addresses the use of swarm intelligence algorithms to achieve a feasible
solution.

##  5 VRFT with robustness constraints

The proposed method regards a two-step procedure. The first step follows the
design of a controller using the VRFT, as commented in Section 3. The same
data from the first step is used in the second step. The need for a second
experiment is avoided since the estimation of $M_{S}$, represented in this
context as $\hat{M}_{S}(\rho)$, does not require new data, as derived in
Subsection 4.1. In this case, the cost function $J^{VR}$ is modified by the
addition of a robustness constraint, regarding the value of the estimated
($\hat{M}_{S}(\rho)$) and desired ($M_{Sd}$) $||S(z,\rho)||_{\infty}$ norm,
leading to a new optimization problem: —l— ρJ^VR(ρ) ^M_S(ρ) ≤M_Sd , which can
be applied directly to the cost function as a penalty [35], resulting in the
Swarm Intelligence optimization cost function: —l— ρJ^SI (ρ) = —— u(k) -
C(z,ρ) (T_d^-1(z) - 1) y(k) ——_2^2 + c H(ρ), where $c$ is a positive constant,
usually with $c\gg 1$, and the $\mathcal{H}_{\infty}$ penalty can be written
as

$H(\rho)=\frac{1}{2}(max[0,\hat{M}_{S}(\rho)-M_{Sd}])^{2}.$ (20)

The controller $C(z,\rho)$ in (5) can be chosen as belonging to any controller
class $\mathcal{C}$, linearly parameterized in $\rho$ as
$\mathcal{C}=\\{C(z,\rho),\rho\in\mathbb{R}^{p}\\}$. The robustness index
$\hat{M}_{S}(\rho)$ is estimated at each iteration of the swarm intelligence
algorithm optimization following the procedure described in Subsection 4.1.

Notice that an important insight can be drawn from (5) regarding the choice of
$M_{Sd}$. If the choice is too ambitious, i.e., $M_{Sd}$ much lower than
$\hat{M}_{S}(\rho_{0})$ estimated with an initial VRFT-obtained controller
with parameter $\rho_{0}$, the obtained parameter $\rho$ could be found too
distant from the VRFT solution. Such a difference could considerably increase
the value of $J^{VR}(\rho)$ to a point where the performance of the swarm
intelligence-obtained controller is drastically affected. On the other hand,
if the constant $c$ is chosen with an exorbitantly high value, only the cases
where the optimization algorithm falls into a local minimum and is not able to
minimize $H(\rho)$ would be affected, resulting in a final cost much higher
than with a lower value of $c$. If the penalty element $H(\rho)$ is solved by
the optimization algorithm, as seen in (20), its value becomes zero,
indifferent of the value of $c$, not affecting the VRFT cost $J^{VR}$.

Considering a search space
$\mathcal{O}\in[l_{b},u_{b}],l_{b},u_{b}\in\mathbb{R}$, in order to accelerate
the convergence of the metaheuristic algorithm, the initialization of the
search agents can inherit the first step solution $\rho_{0}\in\mathbb{R}^{p}$
as a central point, as expressed in

$\overrightarrow{X}_{b}(0)=R\cdot\overrightarrow{X}(0)+\rho_{0},\quad
R=\frac{|\max\\{l_{b},u_{b}\\}|}{2},$ (21)

with $R$ being the initial population spawn radius, and
$\overrightarrow{X}(0)\in\mathbb{R}^{p}$ a random position vector such that
$\overrightarrow{X}(0)\sim U(0,1)$.

An inherent step of the method is to collect input-output data from the
process, as suggested in [27, 11]. Remember to take into account system
identification theory [47] for data to be sufficiently informative. Then, the
two proposed design steps can be applied:

  1. 1.

Use the VRFT to design a controller for the process. Use a flexible reference
model if the plant is NMP, as presented in [34]. Check the obtained robustness
index and proceed to the second step if it does not satisfy $\hat{M}_{S}\leq
M_{Sd}$, or else, take the controller from Step 1 as final. Such a main step
can be divided into the following specific steps:

     * (a)

acquire a data set $\\{u,y\\}_{k=1}^{N}$ from the closed-loop system with an
initial stabilizing controller;

     * (b)

use the data set to design a controller with the VRFT method, as detailed in
Section 3. Controller parameters $\rho$ are obtained after the minimization
procedure of the VRFT method. In the NMP case, parameters for the reference
model $\hat{\eta}$ are also obtained;

     * (c)

estimate the robustness index according to the method in Subsection 4.1. If
$\hat{M}_{S}>M_{Sd}$ proceed to the second step, or else, use the VRFT-
obtained controller with no further modification.

  2. 2.

Apply a swarm intelligence algorithm, as presented next in Subsection 5.1,
considering the optimization problem described in (5) according to the desired
value of $M_{S}$, with restriction applied in the form of a penalty as (20),
with the initial spawn of agents following the recommendation of (21). The
second step can be divided in:

     * (a)

implement the VRFT cost function with the penalty as in (20) regarding the
desired maximum value of $M_{S}$;

     * (b)

change the initialization procedure of the chosen swarm intelligence algorithm
to consider a center spawn $\rho_{0}$, i.e., the VRFT-obtained solution at the
first step, and a spawn radius as suggested in (21) to accelerate convergence;

     * (c)

execute the algorithm and obtain controller parameters that satisfy the
robustness restriction.

The swarm intelligence algorithms considered for obtaining the results of this
work are described next.

###  5.1 Swarm intelligence algorithms

Swarm intelligence algorithms consist of stochastic optimization algorithms
that are based on the collective intelligence of groups composed of simple
agents, usually based on the behavior of animals in nature [57]. In order to
cope with the NFL theorems [38], four algorithms are chosen to be used:
Particle Swarm Optimization (PSO) [43] and Artificial Bee Colony (ABC) [58],
which are well known and widely used [56, 59]; Grey Wolf Optimizer (GWO) [42]
and Improved Grey Wolf Optimizer (I-GWO) [45], more recent approaches with
fewer hyperparameters than the aforementioned. In the following, the four
algorithms are briefly presented.

####  5.1.1 Particle Swarm Optimization

Particle Swarm Optimization (PSO) involves populations (or swarms) in which
each element is called a particle that represents a form of directed mutation
[59, 43]. The swarm is composed of $\ell$ particles searching in a
$D$-dimensional space, which are initialized randomly within the search space
with lower and upper bounds $l_{b}$ and $u_{b}$. Each particle has its own
position ($\overrightarrow{X}_{i}$) and velocity ($\overrightarrow{V}_{i}$)
and is considered as a possible solution for the problem. The best solution
found locally by a particle $i$ is represented by
$\overrightarrow{P}_{i}=\\{P_{i1},P_{i2},...,P_{iD}\\}$, while
$\overrightarrow{G}=\\{G_{1},G_{2},...,G_{D}\\}$ is the best solution found
globally. As for the standard algorithm, each particle is initialized at a
random location with random velocity.

The pseudo-code for PSO is given in A, in which $max\\_it$ is the maximum
number of iterations and $f$ is the cost function. The constant $w_{1}$ is an
inertia weight, $C_{1}$ is the cognitive learning factor, and $C_{2}$ is the
social learning factor. Details for the implementation are available in the
literature [43].

####  5.1.2 Artificial Bee Colony

The Artificial Bee Colony (ABC) algorithm simulates the behavior of bees
performed during their foraging process, conducting a local search in each
iteration. Possible solutions are represented by food sources, while the
quality of each solution is proportional to the nectar amount in each source
[56, 44, 58]. There are three types of bees: scout, employed, and onlooker. At
the initialization, the scout bees randomly find possible food sources
(solutions). Each food source receives an employed bee. By roulette wheel
selection, onlooker bees choose food sources to be exploited based on their
quality, but both types perform local searches in their neighborhood.

The pseudo-code in B describes the ABC optimization procedure, where
$\overrightarrow{X}_{i}$ with $D$ dimensions is the location of food sources.
In this form, $\ell$ is the number of possible solutions (food sources) and
the number of scouts, employed, and onlooker bees are taken as the same as the
number of food sources. $L$ is the abandonment criteria, defined by the
designer. Implementation details are available in [44].

####  5.1.3 Grey Wolf Optimizer

The Grey Wolf Optimizer (GWO) is an algorithm based on the hunting behavior of
grey wolves, which have a strict social dominant hierarchy. The leaders are
the alphas ($\alpha$), who are responsible for making decisions. At the second
level are the betas ($\beta$), subordinates to the alphas who help in
decision-making and other pack activities. The third level wolves are the
deltas ($\delta$), representing scouts, sentinels, elders, hunters, and
caretakers. The rest of the pack is called omega, which must submit to the
higher-ranking wolves [45]. In the GWO algorithm, all wolves follow the mean
position of the $\alpha$, $\beta$, and $\delta$s. At each iteration, the new
three best wolves are re-defined, with the $\alpha$ position being the final
solution. The pseudo-code for the GWO algorithm is presented in C. Details for
the implementation can be checked at [42].

####  5.1.4 Improved Grey Wolf Optimizer

There are three main problems noted in literature around the GWO algorithm
[45]: i) lack of population diversity; ii) imbalance between exploitation and
exploration; iii) premature convergence. The Improved Grey Wolf Optimizer
(I-GWO) tries to solve those issues by changing the search strategy of the GWO
algorithm, including the Dimension Learning-based Hunting (DLH) strategy,
which defines a new (possible) update position for each wolf based not just on
the position of $\alpha$, $\beta$, and $\delta$ wolves, but also on the
position of its neighbors. D briefly describes the I-GWO. Details of the
implementation are given at [45].

##  6 Validation results

In order to validate and illustrate the proposed method, two real-world-
inspired examples are considered. The method is applied as suggested in
Section 5 with all four swarm intelligence algorithms commented in Subsection
5.1. The results are compared in terms of: i) cost value obtained for best
solution - best (lowest) cost; ii) $||S(z,\hat{\rho})||_{\infty}$ value
obtained for best solution; iii) convergence speed. Notice that the system
model is only used to generate data in simulation. The knowledge of the model
is neglected at any stage of the design, maintaining a pure data-driven
fashion.

###  6.1 Example 1: a second-order non-minimum phase plant

The first system to be considered is

$G_{1}(z)=\frac{-0.05(z-1.4)}{z^{2}-1.7z+0.7325}$ (22)

with a time step of 1 second, which is similar to the discrete-time model of a
Boost/Buck-Boost converter operating in continuous conduction mode, regarding
the transfer function of the output voltage by the duty cycle [60]. The
presence of a non-minimum phase zero makes it necessary to use the VRFT with a
flexible criterion [27] at the first step of the proposed method.

Assuming that the system model (22) is unknown, there is no previous knowledge
about its zero being Non-Minimum Phase (NMP). In this sense, it is possible to
analyze the estimated impulse response since the IR of NMP systems initially
moves in the opposite direction (downwards) related to the steady-state one
[51]. Therefore, a Pseudo-Random Binary Signal (PRBS), which is persistently
exciting of high order [47], containing $N=2000$ samples is applied to
$G_{1}(z)$ in simulation, generating an output signal. Additive white Gaussian
noise with a Signal-to-Noise Ratio (SNR) of 20 dB was added to the system at
the output, representing measurement noise. With the input-output dataset, the
IR of $G_{1}(z)$ can be identified with an IR identification algorithm
available in the literature [52, 54, 53, 55], resulting in the signal
presented in Figure 3. Clearly, the IR initially goes downwards, indicating
the presence of an NMP zero, justifying the VRFT with a flexible reference
criterion [27].

Figure 3: Identified impulse response of $G_{1}(z)$.

####  6.1.1 Data collection

The data for estimation is obtained in a closed-loop with a proportional
stabilizing controller [11] since its presence in the system avoids signal
divergence. By the small gain theorem [28], a stabilizing controller can be
obtained following

$k_{p}<\frac{1}{||G||_{\infty}}.$ (23)

Therefore, the stabilizing controller $k_{p}$ is chosen as

$k_{p}=\frac{0.5}{||G_{1}(z)||_{\infty}}=0.8039.$ (24)

In order to obtain the $\mathcal{H}_{\infty}$ norm of $G_{1}(z)$, its impulse
response is estimated according to [54] and the norm is calculated as proposed
in Subsection 4.1.

The excitation signal considered for the data acquisition is a PRBS with
$N=2000$ samples. The signal is applied to the control reference of the
closed-loop system formed by $G_{1}(z)$ with stabilizing controller $k_{p}$.
The control output signal $u(k)$ and the system output signal $y(k)$ are
acquired, forming the input-output set $\\{u,y\\}_{k=1}^{N}$.

####  6.1.2 Step 1 - VRFT with flexible criterion

Assume a situation where the control requirements are: i) null error in
steady-state; ii) settling time approximately 2.5 times faster than in closed-
loop with the stabilizing controller $k_{p}$; iii) null overshoot for a step
reference. A reference model that meets such requirements, following the
guidelines of [11], is chosen as

$T_{d}(z,\hat{\eta_{0}})=\frac{-21(z-1.01)}{(z-0.7)(z-0.3)}.$ (25)

Notice that the initial zero of $T_{d}(z)$ is set as greater than 1, as
suggested in [27], allowing for the VRFT with a flexible criterion to identify
the plant’s NMP zero. The chosen controller class chosen to be used is the PID
class of controllers, which gives

$\bar{C}(z)=\left[1\quad\frac{z}{z-1}\quad\frac{z-1}{z}\right]^{\prime}.$ (26)

After solving the cost function (2) according to the VRFT method with flexible
criterion, the following solution pair for $\eta,\rho$ is obtained:

$\hat{\eta}=\begin{bmatrix}-0.4793&0.6377\end{bmatrix}^{\prime},$ (27a)
$\hat{\rho}=\begin{bmatrix}k_{p}&k_{i}&k_{d}\end{bmatrix}^{\prime}=\begin{bmatrix}1.1246&0.3124&6.9713\end{bmatrix}^{\prime},$
(27b)

resulting in a new reference model $T_{d}(z,\hat{\eta})$, and in the
controller $C(z,\hat{\rho})$, respectively:

$T_{d}(z,\hat{\eta})=\hat{\eta}F(z)=\frac{-0.6899(z-1.33)}{(z-0.7)(z-0.2401)};$
(28a)
$C(z,\hat{\rho})=\hat{\rho}^{\prime}\bar{C}(z)=\frac{8.4083(z^{2}-1.792z+0.8291)}{z(z-1)}.$
(28b)

The non-dominant pole of the reference model, now $T(z,\hat{\eta})$, is
updated together with the minimization of $\eta$ and $\rho$, as suggested in
[11].

By estimating the $\mathcal{H}_{\infty}$ norm of $S(z,\hat{\rho})$ of the
closed-loop system based on the solution (28), an $\hat{M}_{S}=2.1952$ is
obtained, which may be too high for applications that require lower robustness
indexes since it is greater than $2$ [28]. The next subsection presents the
application of the second step of the proposed method to reduce $M_{S}$ for
the obtained VRFT solution.

####  6.1.3 Step 2 - Swarm intelligence algorithm

The use of swarm intelligence algorithms for solving the proposed problem is
straightforward. At their implementation - see A, B, C, and D for the pseudo-
codes - the cost function $f$ is the proposed $J^{SI}(\rho)$ (5). At each
iteration, the input to $f$ is the parameter $\rho$, obtained in the previous
iteration for each search agent since it is required to parameterize the
controller $C(z,\rho)$ and to estimate the $\mathcal{H}_{\infty}$ norm
$\hat{M}_{S}(\rho)$.

The computation of $J^{SI}(\rho)$ is done by directly calculating the cost
(5). The norm $\hat{M}_{S}(\rho)$ can be obtained through (15), estimating the
impulse response of $S(z,\rho)$ with the signals proposed in (19) and an IR
estimation algorithm available in the literature, such as [54].

The upper search bound for the swarm intelligence algorithms is defined as
$u_{b}=10$, which should be sufficient considering that the maximum $\rho$
value of the VRFT-obtained controller is $6.9713$ and, taking into account a
choice of $M_{Sd}$ that is not too ambitious, the resultant $\rho$ should not
be too distant from the initialized value. The lower search bound is chosen as
$l_{b}=0$ to avoid negative controller gain, making the obtained controller
passive [61].

The initialization of all agents is done randomly with a uniform distribution
within the suggested spawn radius in (21), $R=u_{b}/2=5$, with the central
point equal to the VRFT solution at Step 1 (27b). The reference model for the
cost function (5) is considered to be the VRFT with a flexible criterion
reference model, $T_{d}(z,\hat{\eta)}$, obtained in Step 1 as (28a). The
desired $M_{S}$ to be achieved, $M_{Sd}$, is set to 1.8, which is a sufficient
value in terms of robustness, satisfies $M_{Sd}\leq 2$ and should not
compromise substantially the performance of the system, which could happen if
$M_{Sd}<<2$.

To make the comparison between algorithms possible, the number of agents was
fixed to 50 and the number of iterations was limited to 100. In order to
obtain a satisfactory number of realizations for the analysis of results each
algorithm was executed 50 times. The hyperparameters for the PSO and ABC
algorithms, except for the number of agents and the maximum number of
iterations, are presented in Table 1. The PSO parameters were chosen as the
MATLAB® default parameters of the Global Optimization Toolbox [53], while the
ABC parameters were used according to the algorithm implementation of [62].
GWO and I-GWO do not contain any hyperparameter aside from the number of
agents and the maximum number of iterations.

Algorithm | Parameter settings | Value
---|---|---
PSO | Cognitive learning factor ($C_{1}$) | 1.49
Social learning factor ($C_{2}$) | 1.49
Inertia range (range of $w_{1}$) | [0.1,1.1]
ABC | Abandonment criteria ($L$) | 90
Acceleration coefficient ($a$) | 1
Table 1: Parameters settings for PSO and ABC.

Figure 4 shows the average convergence curve of all algorithms for 50 runs,
considering the system $G_{1}(z)$ as aforementioned. Table 2 presents the time
that each iteration took and the number of iterations to converge, considering
the average value for all 50 realizations and a convergence criterion of
$\delta=1\times 10^{-3}$ from an iteration cost to the subsequent one. The
results were obtained with an Intel Core i5 4670 3.40 GHz processor. The I-GWO
algorithm took a longer time to converge, followed by ABC, PSO, and finally,
GWO. Comparing similar algorithms, I-GWO obtained a duration for one iteration
that is more than twice the same duration for the GWO. But since the
optimization of the proposed cost function is executed offline, where the
obtained parameters ($\rho$) are applied to the controller without any further
modification, this is not considered an issue for the current application.

Considering all 50 executions per swarm intelligence algorithm, Figure 5
presents the best (lower) cost statistics for all algorithms at example 1 in
the form of a box plot. Clearly, I-GWO achieved the most desirable performance
in terms of cost since it contains fewer outliers and a very low dispersion if
compared to the other algorithms’ solutions. PSO, ABC, and GWO, in general,
result in higher cost values than I-GWO for the considered cost function.
Table 3 shows the quantitative values related to the best cost of all
algorithms at each run, confirming the conclusions taken from Figure 5. The
obtained values for the parameter $\rho=[k_{p}\ k_{i}\ k_{d}]^{\prime}$ are
presented in Figure 6 and show how the lower cost dispersion obtained with the
I-GWO algorithm can influence the values of the controller parameters and
their tendency. For the case of example 1, the I-GWO is the most robust
algorithm in relation to the random initialization of the parameters. Such a
characteristic is desired in a real-world scenario, reducing the number of
trials in the design procedure to obtain a satisfactory solution. Also, notice
that the median solution obtained with the I-GWO is not too far from the VRFT
solution (27b) obtained at Step 1, which agrees with the previous hypothesis
that a value of desired robustness $M_{Sd}$ chosen as not too ambitious when
compared with the initial VRFT solution (27b) should result in a final
solution that is close to the initial (VRFT) one.

Figure 4: Average convergence curves for all algorithms considering a Monte Carlo experiment of 50 executions for example 1. Algorithm | 1-it. time (s) | It. to converge | Time to converge (s)
---|---|---|---
PSO | $\ 8.04$ | $41$ | $329.75$
ABC | $22.63$ | $19$ | $430.02$
GWO | $11.25$ | $14$ | $157.53$
I-GWO | $24.18$ | $20$ | $483.66$

Table 2: Time for convergence of all algorithms for example 1.
Figure 5: Box plot of a Monte Carlo experiment with 50 executions for all algorithms in terms of best cost value obtained for example 1. Algorithm | median | $\sigma$ | min | max
---|---|---|---|---
PSO | $0.2032$ | $0.2782$ | $0.2017$ | $2.1710$
ABC | $0.2334$ | $0.0402$ | $0.2025$ | $0.4009$
GWO | $0.2489$ | $0.0858$ | $0.2017$ | $0.4764$
I-GWO | $0.2017$ | $5.0516\times 10^{-5}$ | $0.2017$ | $0.2018$
Table 3: Quantitative results from the box plot in terms of best cost for
example 1.

(a) $k_{p}$

(b) $k_{i}$

(c) $k_{d}$

Figure 6: Box plot of the obtained $k_{p}$ ($\rho_{1}$), $k_{i}$ ($\rho_{2}$),
and $k_{d}$ ($\rho_{3}$) controller parameters for all algorithms evaluated at
50 executions for example 1.

Figure 7 presents the box plot for the obtained $||S(z,\hat{\rho})||_{\infty}$
by the best solution of each algorithm at each run, in a closed-loop with
$G_{1}(z)$. I-GWO obtained the most desired result in terms of $\hat{M}_{S}$
considering the lack of outliers and the low dispersion. PSO had one outlier
with $\hat{M}_{S}>2$, while ABC obtained three outliers of higher
$\hat{M}_{S}$, and the performance by the GWO algorithm for this problem was
not satisfactory since there were too many solutions that achieved an
$\hat{M}_{S}$ higher than 2. Table 4 shows the quantitative data of the box
plot presented in Figure 7, in agreement with what is commented over the
results.

Figure 7: Box plot of 50 executions for all algorithms in terms of $||S(z,\hat{\rho})||_{\infty}$ value obtained for example 1. Algorithm | median | $\sigma$ | min | max
---|---|---|---|---
PSO | $1.8014$ | $0.1130$ | $1.7509$ | $2.5988$
ABC | $1.8091$ | $0.4740$ | $9.9798\times 10^{-5}$ | $2.8496$
GWO | $1.8087$ | $0.2923$ | $1.7995$ | $2.5760$
I-GWO | $1.8020$ | $4.5561\times 10^{-4}$ | $1.8008$ | $1.8030$
Table 4: Quantitative results from the box plot in terms of
$||S(z,\hat{\rho})||_{\infty}$ for example 1.

To briefly demonstrate that, although the robustness is increased with the
proposed method, the reference tracking is not overly penalized, a step
reference was applied to the controlled plant of example 1 using a controller
designed by the proposed method with the I-GWO algorithm, with parameters
$\rho_{IGWO}=[1.090\quad 0.2194\quad 5.4018]$, and the VRFT controller (27b)
obtained at step 1. As shown in Figure 8, the VRFT-controlled system achieves
a settling time of 54 seconds, a step overshoot of 20 % and an undershoot of
42 %. For the proposed method, a settling time of 39 seconds is achieved, with
9 % of overshoot and 33 % of undershoot. Notice that, in this specific case,
the proposed method could even enhance the reference tracking performance in
terms of settling time, however, this is not expected in most of the systems
since the inclusion of a robustness constraint will penalize the VRFT cost
function. The reduced overshoot and undershoot are a consequence of the
increased robustness, with an $\hat{M}_{S}$ of 2.2079 for the VRFT and 1.8030
for the I-GWO-designed controller.

Figure 8: Output signal for a reference step signal applied to the controlled
plant with a controller designed with the proposed method using I-GWO and the
VRFT-designed controller for example 1.

###  6.2 Example 2: fourth-order plant

The fourth-order plant consists of

$G_{2}(z)=\frac{0.1381(z-0.95)(z^{2}-1.62z+0.6586)}{(z^{2}-1.7z+0.7325)(z^{2}-1.84z+0.8564)},$
(29)

with a time step of 1 second, which has the same structure as the discrete-
time transfer function of a SEPIC converter, from duty cycle to output voltage
[63]. Since the plant’s zeros have minimum phase, which can be evaluated with
data as aforementioned in Subsection 6.1, the VRFT method is used without the
flexible reference model [27].

####  6.2.1 Data collection

For plant $G_{2}(z)$, the data is obtained in the same way as described, for
example 1, in Subsection 6.1, with a PRBS signal of $N=2000$ samples applied
to the closed-loop system with stabilizing controller

$k_{p}=\frac{0.5}{||G_{2}(z)||_{\infty}}=0.3828,$ (30)

considering additive white Gaussian noise with an SNR of 20 dB to represent
measurement noise. The input-output set is formed by $\\{u,y\\}_{k=1}^{N}$.

####  6.2.2 Step 1 - VRFT

After the data is acquired, the next step is to use the VRFT to design a
controller, which solves the cost function (2). For this example, the
following control requirements are assumed: i) null steady-state error; ii)
settling time of approximately 6.5 times faster than the closed-loop settling
time with stabilizing controller $k_{p}$; iii) null overshoot for a step
reference. Considering such requirements, the choice of the reference model is
made as suggested in [11], obtaining

$Td(z)=\frac{1.4(z-0.6)}{(z-0.3)(z-0.2)}.$ (31)

Supposing a limited situation where only a PI controller is available, e.g.,
because of hardware limitations on a certain product. Therefore, the
controller class to be considered is the PI class of controllers, resulting in

$\bar{C}(z)=\left[1\quad\frac{z}{z-1}\right]^{\prime}.$ (32)

The obtained VRFT solution results in the controller parameter

$\hat{\rho}=[k_{p}\quad k_{i}]=[6.6568\quad 3.3728].$ (33)

Via (3), the controller

$C(z,\hat{\rho})=\hat{\rho}^{\prime}\bar{C}(z)=\frac{10.03(z-0.6637)}{(z-1)}$
(34)

is obtained.

Considering the VRFT-obtained solution (33), the robustness index of the
system can be estimated according to Subsection 4, obtaining
$\hat{M}_{S}=2.2767$. As aforementioned, an $M_{S}\leq 2$ is desired to ensure
sufficient robustness [28], which leads to the application of the second step
of the proposed method.

####  6.2.3 Step 2 - Swarm intelligence algorithm

The swarm intelligence algorithms PSO, ABC, GWO, and I-GWO are applied to the
problem (5) for the fourth-order plant case, according to the pseudo-codes
presented in A, B, C, and D, where $f$ is the cost function (5), as already
commented at Subsection 6.1.3. The upper search bound is kept at $u_{b}=10$
and the lower bound at $l_{b}=0$, in order to increase the passivity of the
controller as mentioned in Subsection 6.1. An upper bound of $10$ should be
sufficient, considering that the maximum desired robustness is not too far
from the estimated robustness index at the end of step 1. The initial
population spawn radius follows (21), $R=(|u_{b}|+|l_{b}|)/2=5$, with the
solution of the VRFT at Step 1 (33) as the central point. The desired
$||S(z,\hat{\rho})||_{\infty}$ is set to 1.5, which satisfies $M_{Sd}\leq 2$
and is not much lower than 2.

The number of agents of all algorithms is set to 50, with a maximum of 100
iterations per execution. Each algorithm is executed 50 times for different
noise realizations, such that the results are sufficiently representative for
further analysis. For PSO and ABC algorithms, parameters are set as presented
in Table 1. Aside from the number of agents and the maximum number of
iterations, no other parameter is set by the user with the proposed GWO and
I-GWO algorithm.

The average convergence curve of all algorithms for this case is presented in
Figure 9. Table 5 presents the time for one iteration and the number of
iterations each algorithm took to converge, considering a convergence
criterion $\delta=1\times 10^{-3}$. The hardware configuration is the same as
in example 1. The ABC algorithm was the slowest algorithm in this example,
followed by I-GWO, PSO, and at last, GWO. In comparison to example 1, GWO is
still the fastest algorithm to converge, while ABC took I-GWO place as the
slowest convergence.

Figure 10 shows the box plot regarding the best fitness value for each
algorithm, considering all executions. PSO, ABC, and I-GWO did not present far
outliers, as those seen in the case of GWO. The quantitative values of the box
plot are shown in Table 6. In this case, the most desirable result is obtained
by the PSO algorithm. Although, notice that the scale of the y-axis (best
cost) only varies at the third decimal place, which means that PSO, ABC, and
I-GWO solutions should perform very similarly in practice.

Figure 9: Average convergence curves for all algorithms considering a Monte Carlo experiment of 50 executions for example 2. Alg. | 1-it. time (s) | It. to converge | Time to converge (s)
---|---|---|---
PSO | $7.24$ | $20$ | $144.72$
ABC | $22.93$ | $10$ | $229.32$
GWO | $7.55$ | $9$ | $68.00$
I-GWO | $22.34$ | $9$ | $201.10$

Table 5: Time for convergence of all algorithms for example 2.
Figure 10: Box plot of a Monte Carlo experiment with 50 executions for all algorithms in terms of best cost value obtained for example 2. Algorithm | median | $\sigma$ | min | max
---|---|---|---|---
PSO | $0.49284$ | $2.2270\times 10^{-9}$ | $0.49284$ | $0.49284$
ABC | $0.49287$ | $4.7174\times 10^{-5}$ | $0.49284$ | $0.49303$
GWO | $0.49290$ | $1.0414\times 10^{-2}$ | $0.49284$ | $0.53568$
I-GWO | $0.49285$ | $6.4636\times 10^{-6}$ | $0.49284$ | $0.49286$
Table 6: Quantitative results from the box plot in terms of best cost for the
example with system $G_{2}(z)$.

The $||S(z,\hat{\rho})||_{\infty}$ norm obtained for the best solution at each
run is shown in Figure 11, with its quantitative values presented in Table 7.
PSO, ABC, and I-GWO present very similar values for the obtained
$\mathcal{H}_{\infty}$ norm, which shows that the performance of those
solutions is close, as also noticed from the cost analysis. GWO is the only
algorithm that presents outliers in this case, which, in comparison to I-GWO,
indicates that the local minima problems pointed out in literature [45] indeed
can occur.

Additionally, similar conclusions can be drawn from the controller parameter’s
values obtained with each optimization algorithm for the proposed problem. In
Figure 12, outliers are only obtained with the GWO algorithm, where the gain
of the integral part of the controller, $k_{i}$, achieves values close to
zero, which has, as a consequence, a lower performance for the controlled
system, even though the robustness criteria is met. The obtained $\rho$ for
PSO, ABC, and I-GWO are very similar for all executions, indicating that any
of the three algorithms could be used for this example without significant
loss of performance or robustness when compared to one another.

Figure 11: Box plot of a Monte Carlo experiment with 50 executions for all algorithms in terms of $||S(z,\hat{\rho})||_{\infty}$ for example 2. Algorithm | median | $\sigma$ | min | max
---|---|---|---|---
PSO | $1.5169$ | $1.4558\times 10^{-6}$ | $1.5169$ | $1.5169$
ABC | $1.5165$ | $4.4060\times 10^{-4}$ | $1.5148$ | $1.5172$
GWO | $1.5133$ | $8.4447\times 10^{-3}$ | $1.4825$ | $1.5173$
I-GWO | $1.5168$ | $1.2216\times 10^{-4}$ | $1.5165$ | $1.5171$
Table 7: Quantitative results from the box plot in terms of
$||S(z,\hat{\rho})||_{\infty}$ for example 2.

(a) $k_{p}$

(b) $k_{i}$

Figure 12: Box plot of the obtained $k_{p}$ ($\rho_{1}$) and $k_{i}$
($\rho_{2}$) controller parameters for all algorithms evaluated at 50
executions for example 2.

In the same way as for example 1, a step reference is applied to the
controlled plant of example 2 using a VRFT-designed controller and a
controller designed with the proposed method using I-GWO, as shown in Figure
13. In this case, as one would expect, the reference tracking was penalized in
terms of settling time, from 5 seconds (VRFT) to 14 seconds (I-GWO-designed).
The overshoot is 42 % for the VRFT controller and 29 % for the proposed method
with I-GWO, reflecting the reduced $\hat{M}_{S}$ \- from 2.3834 to 1.5156.

Figure 13: Output signal for a reference step signal applied to the controlled
plant with a controller designed with the proposed method using I-GWO and the
VRFT-designed controller for example 2.

##  7 Conclusion

This work proposed a data-driven one-shot technique to increase the robustness
of a closed-loop discrete-time system by changing the controller parameters
using swarm intelligence algorithms. The considered optimization problem (5)
is the VRFT cost function with the addition of a penalty regarding the value
of the $||S(z,\rho)||_{\infty}$ norm, which can be directly used as a measure
of robustness. Such a value is estimated via an impulse response at each
iteration of the metaheuristic algorithm.

Four swarm intelligence algorithms - PSO, ABC, GWO, and I-GWO - have been
considered to illustrate the proposed technique with two real-world inspired
plants. For the example of a second-order non-minimum phase plant, I-GWO
obtained the best results in terms of dispersion, outliers, and cost, with
acceptable values of $||S(z,\rho)||_{\infty}$. PSO achieved similar results in
this case, while ABC and GWO algorithms had a higher occurrence of outliers -
in terms of cost and $\hat{M}_{S}$. In the case of the fourth-order minimum-
phase plant, ABC, I-GWO, and PSO obtained very similar results, all
satisfactory in terms of cost and $\mathcal{H}_{\infty}$ norm of $S(z,\rho)$.
GWO, on the other hand, obtained a few outliers that, although the criterion
of robustness was met, those cases achieved a very low integral gain, which
could drastically affect the performance of the system. In terms of speed of
convergence, I-GWO and ABC were the slowest algorithms, while PSO converged
faster and GWO was the fastest.

As for future works, it is suggested: the inclusion of other constraints
(e.g., for control effort) simultaneously with the robustness constraints; the
use of other types of metaheuristics, as evolutionary of physics-based
algorithms; the inclusion of a robustness constraint to the OCI, VDFT, or DD-
LQR methods; extension of the current work for MIMO systems.

##  8 Acknowledgments

This study was financed in part by the Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001, partly by the
Fundação de Amparo à Pesquisa e Inovação do Estado de Santa Catarina (FAPESC)
- Grant number 288/2021, and partly by Conselho Nacional de Desenvolvimento
Científico e Tecnológico - CNPq - Brazil.

## Appendix A

input : $f$, $l_{b}$, $u_{b}$, $\ell$, $max\\_it$, $w_{1}$, $C_{1}$, $C_{2}$

output : $\overrightarrow{G}$

1 initialize $\ell$ particles with random position $\overrightarrow{X}_{i}(0)$
within $l_{b}$, $u_{b}$, and velocity $\overrightarrow{V}_{i}(0)$;

2 define (initial) best local solution $\overrightarrow{P}_{i}$ as
$f[\overrightarrow{X}_{i}(0)]$ for each particle;

3 define (initial) best global solution $\overrightarrow{G}$ as the position
of the particle with the best fitness;

4 for  _$n\in\\{1,...,max\\_it\\}$_ do

5 for  _$i\in\\{1,...,\ell\\}$_ do

6 define new velocities $\overrightarrow{V}_{i}(n)$ for each particle with
$w_{1}$, $C_{1}$, $C_{2}$;

7 define new positions $\overrightarrow{X}_{i}(n)$ for each particle;

8

9 update best local solution for all particles;

10 $\overrightarrow{G}\leftarrow$ position of the particle with higher fitness

Algorithm 1 Particle swarm optimization pseudo-code

## Appendix B

input : $f$, $l_{b}$, $u_{b}$, $\ell$, $max\\_it$, $L$

output : food source with more nectar (higher fitness)

11 initialize $l$ food sources with position $\overrightarrow{X}_{i}(0)$
within boundaries $l_{b}$, $u_{b}$;

12 define (initial) best solution $\overrightarrow{G}$ as the food source with
higher fitness;

13 for  _$n\in\\{1,...,max\\_it\\}$_ do

14 for  _$i\in\\{1,...,\ell\\}$)_ do

15 define a new position for the $i$-th employed bee based on the current food
sources;

16 if  _the new position has higher fitness than the current_ then

17 set as the new position for the employed bee;

18

19 else

20 increase abandonment counter;

21

22

23 for  _$i\in\\{1,...,\ell\\}$_ do

24 select a food source by roulette wheel selection based on the probability
of such position to have more nectar in comparison with others;

25 define a position for the $i$-th onlooker bee based on the selected food
source;

26 if  _the new position has higher fitness than the current_ then

27 set as the new position for the employed bee;

28

29 else

30 increase abandonment counter;

31

32

33 if  _abandonment counter exceeds $L$_ then

34 the exceeding food sources are abandoned and new (random) food sources are
initialized by scout bees;

35 the abandonment counter is reset;

36

37

Algorithm 2 Artificial bee colony

## Appendix C

input : $f$, $l_{b}$, $u_{b}$, $\ell$, $max\\_it$

output : position of the $\alpha$ wolf

38 initialize $\ell$ wolves with random position within boundaries $l_{b}$,
$u_{b}$;

39 define the three higher fitness wolves as $\alpha$, $\beta$, and $\delta$;

40 for  _$n\in\\{1,...,max\\_it\\}$_ do

41 for  _$i\in\\{1,...,\ell\\}$_ do

42 position is updated following the mean of the $\alpha$, $\beta$, and
$\delta$ wolves;

43

44 the three wolves with higher fitness are the new $\alpha$, $\beta$, and
$\delta$;

45

Algorithm 3 Grey wolf optimizer pseudo-code

## Appendix D

input : $f$, $l_{b}$, $u_{b}$, $\ell$, $max\\_it$

output : position of the $\alpha$ wolf

46 initialize $\ell$ wolves with random position within boundaries $l_{b}$,
$u_{b}$;

47 define the three higher fitness wolves as $\alpha$, $\beta$, and $\delta$;

48 for  _$n\in\\{1,...,max\\_it\\}$_ do

49 for  _$i\in\\{1,...,\ell\\}$_ do

50 the candidate position is defined following the mean of the $\alpha$,
$\beta$, and $\delta$ wolves (as in GWO);

51 a neighborhood is constructed based on the current position and candidate
position;

52 multi-neighbor learning is performed (DLH solution);

53 choose best position between candidate and DLH for the update;

54

55 the three wolves with higher fitness are the new $\alpha$, $\beta$, and
$\delta$;

56

Algorithm 4 Improved grey wolf optimizer pseudo-code

## References

  * [1] N. S. Nise, Control Systems Engineering, 3rd Edition, John Wiley & Sons, Inc., USA, 2000.
  * [2] N. Chaudhuri, D. Chakraborty, B. Chaudhuri, Damping control in power systems under constrained communication bandwidth: A predictor corrector strategy, IEEE Transactions on Control Systems Technology 20 (1) (2012) 223–231, cited By 21.  doi:10.1109/TCST.2010.2096817.
  * [3] R. Xie, I. Kamwa, C. Y. Chung, A novel wide-area control strategy for damping of critical frequency oscillations via modulation of active power injections, IEEE Transactions on Power Systems 36 (1) (2021) 485–494.  doi:10.1109/TPWRS.2020.3006438.
  * [4] M. Kazimierczuk, Pulse-width Modulated DC-DC Power Converters, Wiley, 2008.
URL https://books.google.com.br/books?id=H_FwBkSpYY0C

  * [5] T. Kobaku, S. C. Patwardhan, V. Agarwal, Experimental evaluation of internal model control scheme on a dc–dc boost converter exhibiting nonminimum phase behavior, IEEE Transactions on Power Electronics 32 (11) (2017) 8880–8891.  doi:10.1109/TPEL.2017.2648888.
  * [6] R. A. Aguiar, I. C. Franco, F. Leonardi, F. Lima, Fractional pid controller applied to a chemical plant with level and ph control, Chemical Product and Process Modeling 13 (2018).
  * [7] V. Tharanidharan, R. Sakthivel, Y. Ren, S. Marshal Anthoni, Robust finite-time pid control for discrete-time large-scale interconnected uncertain system with discrete-delay, Mathematics and Computers in Simulation 192 (2022) 370–383.  doi:10.1016/j.matcom.2021.08.024.
  * [8] J. Tudon-Martinez, J.-J. Lozoya-Santos, A. Cantu-Perez, A. Cardenas-Romero, Advanced temperature control applied on an industrial box furnace, Journal of Thermal Science and Engineering Applications 14 (6), cited By 0 (2022).  doi:10.1115/1.4052020.
  * [9] N. van Tan, K. Dang, P. Dai, L. Van, Position control for haptic device based on discrete-time proportional integral derivative controller, International Journal of Electrical and Computer Engineering 12 (1) (2022) 269–276, cited By 0.  doi:10.11591/ijece.v12i1.pp269-276.
  * [10] L. H. Keel, S. P. Bhattacharyya, Controller synthesis free of analytical models: Three term controllers, IEEE Transactions on Automatic Control 53 (6) (2008) 1353–1369.  doi:10.1109/TAC.2008.925810.
  * [11] C. L. Remes, R. B. Gomes, J. V. Flores, F. B. Líbano, L. Campestrini, Virtual reference feedback tuning applied to dc–dc converters, IEEE Transactions on Industrial Electronics 68 (1) (2021) 544–552.  doi:10.1109/TIE.2020.2967729.
  * [12] I. Zenelis, X. Wang, A model-free sparse wide-area damping controller for inter-area oscillations, International Journal of Electrical Power and Energy Systems 136, cited By 0 (2022).  doi:10.1016/j.ijepes.2021.107609.
  * [13] H. Huang, X. Huang, W. Ding, M. Yang, D. Fan, J. Pang, Uncertainty optimization of pure electric vehicle interior tire/road noise comfort based on data-driven, Mechanical Systems and Signal Processing 165, cited By 0 (2022).  doi:10.1016/j.ymssp.2021.108300.
  * [14] H. Hjalmarsson, M. Gevers, S. Gunnarsson, O. Lequin, Iterative feedback tuning: theory and applications, IEEE Control Systems Magazine 18 (4) (1998) 26–41.  doi:10.1109/37.710876.
  * [15] A. Karimi, L. Mišković, D. Bonvin, Iterative correlation-based controller tuning, International Journal of Adaptive Control and Signal Processing 18 (8) (2004) 645–664.  doi:https://doi.org/10.1002/acs.825.
  * [16] J. Li, Y. Fang, Y. Ge, Y. Wu, Dnn-based implementation of data-driven iterative learning control for unknown system dynamics, in: 2020 IEEE 9th Data Driven Control and Learning Systems Conference (DDCLS), 2020, pp. 1037–1042.  doi:10.1109/DDCLS49620.2020.9275089.
  * [17] D. Meng, J. Zhang, Design and analysis of data-driven learning control: An optimization-based approach, IEEE Transactions on Neural Networks and Learning Systems 33 (10) (2022) 5527–5541.  doi:10.1109/TNNLS.2021.3070920.
  * [18] M. Campi, A. Lecchini, S. Savaresi, Virtual reference feedback tuning: a direct method for the design of feedback controllers, Automatica 38 (8) (2002) 1337–1346.  doi:https://doi.org/10.1016/S0005-1098(02)00032-8.
URL https://www.sciencedirect.com/science/article/pii/S0005109802000328

  * [19] L. C. Kammer, R. R. Bitmead, P. L. Bartlett, Direct iterative tuning via spectral analysis, Automatica 36 (9) (2000) 1301–1307.  doi:https://doi.org/10.1016/S0005-1098(00)00040-6.
  * [20] L. Campestrini, D. Eckhard, A. Sanfelice Bazanella, M. Gevers, Data-driven model reference control design by prediction error identification, Journal of the Franklin Institute 354 (6) (2017) 2628–2647, special issue on recent advances on control and diagnosis via process measurements.  doi:https://doi.org/10.1016/j.jfranklin.2016.08.006.
  * [21] D. Eckhard, L. Campestrini, E. C. Boeira, Virtual disturbance feedback tuning, IFAC Journal of Systems and Control 3 (2018) 23–29, n/a.  doi:https://doi.org/10.1016/j.ifacsc.2018.01.003.
  * [22] G. R. GONÇALVES DA SILVA, A. S. Bazanella, C. Lorenzini, L. Campestrini, Data-driven lqr control design, IEEE Control Systems Letters 3 (1) (2019) 180–185.  doi:10.1109/LCSYS.2018.2868183.
  * [23] P. G. Carlet, A. Favato, S. Bolognani, F. Dörfler, Data-driven predictive current control for synchronous motor drives, in: 2020 IEEE Energy Conversion Congress and Exposition (ECCE), 2020, pp. 5148–5154.  doi:10.1109/ECCE44975.2020.9235958.
  * [24] K. Prag, M. Woolway, T. Celik, Toward data-driven optimal control: A systematic review of the landscape, IEEE Access 10 (2022) 32190–32212.  doi:10.1109/ACCESS.2022.3160709.
  * [25] J. A. R. Pérez, R. S. Llopis, Tuning and robustness analysis of event-based pid controllers under different event-generation strategies, International Journal of Control 91 (7) (2018) 1567–1587.  doi:10.1080/00207179.2017.1322716.
URL https://doi.org/10.1080/00207179.2017.1322716

  * [26] S. Alcántara, R. Vilanova, C. Pedret, Pid control in terms of robustness/performance and servo/regulator trade-offs: A unifying approach to balanced autotuning, Journal of Process Control 23 (4) (2013) 527–542.  doi:https://doi.org/10.1016/j.jprocont.2013.01.003.
URL https://www.sciencedirect.com/science/article/pii/S0959152413000139

  * [27] A. S. Bazanella, L. Campestrini, D. Eckhard, Data-Driven Controller Design: The H2 Approach, Springer Publishing Company, Incorporated, 2014.
  * [28] S. Skogestad, I. Postlethwaite, Multivariable feedback control: Analysis and Design, John Wiley, Hoboken, US-NJ, 2005.
  * [29] S. Chiluka, S. Ambati, M. Seepana, U. Babu Gara, A novel robust virtual reference feedback tuning approach for minimum and non-minimum phase systems, ISA Transactions 115 (2021) 163–191.  doi:10.1016/j.isatra.2021.01.018.
  * [30] A. Karimi, C. Kammer, A data-driven approach to robust control of multivariable systems by convex optimization, Automatica 85 (2017) 227–233.  doi:https://doi.org/10.1016/j.automatica.2017.07.063.
  * [31] A. Nicoletti, M. Martino, A. Karimi, A robust data-driven controller design methodology with applications to particle accelerator power converters, IEEE Transactions on Control Systems Technology 27 (2) (2019) 814–821.  doi:10.1109/TCST.2017.2783346.
  * [32] J. Na, J. Zhao, G. Gao, Z. Li, Output-feedback robust control of uncertain systems via online data-driven learning, IEEE Transactions on Neural Networks and Learning Systems 32 (6) (2021) 2650–2662.  doi:10.1109/TNNLS.2020.3007414.
  * [33] J. Berberich, J. Köhler, M. A. Müller, F. Allgöwer, Data-driven model predictive control with stability and robustness guarantees, IEEE Transactions on Automatic Control 66 (4) (2021) 1702–1717.  doi:10.1109/TAC.2020.3000182.
  * [34] L. Campestrini, D. Eckhard, M. Gevers, A. Bazanella, Virtual reference feedback tuning for non-minimum phase plants, Automatica 47 (8) (2011) 1778–1784.  doi:https://doi.org/10.1016/j.automatica.2011.04.002.
  * [35] D. G. Luenberger, Y. Ye, Linear and Nonlinear Programming, Springer Publishing Company, Incorporated, 2015.
  * [36] L. Jiaqi, L. Shuangquan, L. Ning, G. Xiaozhi, Parametric data-driven control method of hybrid system via virtual reference feedback tuning, in: 2017 36th Chinese Control Conference (CCC), 2017, pp. 9676–9681.  doi:10.23919/ChiCC.2017.8028899.
  * [37] D. Selvi, D. Piga, A. Bemporad, Towards direct data-driven model-free design of optimal controllers, in: 2018 European Control Conference (ECC), 2018, pp. 2836–2841.  doi:10.23919/ECC.2018.8550184.
  * [38] D. Wolpert, W. Macready, No free lunch theorems for optimization, IEEE Transactions on Evolutionary Computation 1 (1) (1997) 67–82.  doi:10.1109/4235.585893.
  * [39] S. Mirjalili, Genetic Algorithm, Springer International Publishing, Cham, 2019.
  * [40] B. Alatas, U. Can, Physics based metaheuristic optimization algorithms for global optimization, American Journal of Information Science and Computer Engineering (01 2015).
  * [41] M. N. Ab Wahab, S. Nefti-Meziani, A. Atyabi, A comprehensive review of swarm optimization algorithms, PLOS ONE 10 (05 2015).  doi:10.1371/journal.pone.0122827.
URL https://doi.org/10.1371/journal.pone.0122827

  * [42] S. Mirjalili, S. M. Mirjalili, A. Lewis, Grey wolf optimizer, Advances in Engineering Software 69 (2014) 46–61.  doi:https://doi.org/10.1016/j.advengsoft.2013.12.007.
  * [43] J. Kennedy, R. Eberhart, Particle swarm optimization, in: Proceedings of ICNN’95 - International Conference on Neural Networks, Vol. 4, 1995, pp. 1942–1948 vol.4.  doi:10.1109/ICNN.1995.488968.
  * [44] D. Karaboga, An idea based on honey bee swarm for numerical optimization, Tech. rep., Erciyes University, Engineering Faculty (2005).
  * [45] M. H. Nadimi-Shahraki, S. Taghian, S. Mirjalili, An improved grey wolf optimizer for solving engineering problems, Expert Systems with Applications 166 (2021) 113917.  doi:https://doi.org/10.1016/j.eswa.2020.113917.
  * [46] G. C. Goodwin, S. F. Graebe, M. E. Salgado, Control System Design, 1st Edition, Prentice Hall PTR, USA, 2000.
  * [47] L. Ljung, System Identification (2nd Ed.): Theory for the User, Prentice Hall PTR, USA, 1999.
  * [48] L. V. Fiorio, Virtual reference feedback tuning with robustness constraints: a swarm intelligence solution, Master’s thesis, Universidade do Estado de Santa Catarina (2022).  doi:10.13140/RG.2.2.19418.41927.
  * [49] G. R. Gonçalves da Silva, A. S. Bazanella, L. Campestrini, One-shot data-driven controller certification, ISA transactions 99 (2020) 361—373.  doi:10.1016/j.isatra.2019.10.011.
  * [50] T. Chen, H. Ohlsson, L. Ljung, On the estimation of transfer functions, regularizations and gaussian processes—revisited, Automatica 48 (8) (2012) 1525–1535.  doi:https://doi.org/10.1016/j.automatica.2012.05.026.
  * [51] S. L. Brunton, J. N. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control, 1st Edition, Cambridge University Press, USA, 2019.
  * [52] T. Chen, L. Ljung, Implementation of algorithms for tuning parameters in regularized least squares problems in system identification, Automatica 49 (7) (2013) 2213–2220.  doi:https://doi.org/10.1016/j.automatica.2013.03.030.
  * [53] The MathWorks Inc., Natick, Massachussetts, United States, System Identification Toolbox Release 2017b (2021).
  * [54] L. V. Fiorio, C. L. Remes, Y. R. de Novaes, impulseest: A python package for non-parametric impulse response estimation with input–output data, SoftwareX 15 (2021) 100761.  doi:https://doi.org/10.1016/j.softx.2021.100761.
  * [55] S. Yerramilli, A. Tangirala, sysid: system identification in R, R Foundation for Statistical Computing (2017).
  * [56] K.-L. Du, M. N. S. Swamy, Search and Optimization by Metaheuristics: Techniques and Algorithms Inspired by Nature, 1st Edition, Birkhäuser Basel, 2016.
  * [57] E. Bonabeau, M. Dorigo, G. Theraulaz, Swarm Intelligence: From Natural to Artificial Systems, Oxford University Press, Santa Fe Institute Studies in the Sciences of Complexity, 2001.
  * [58] D. Karaboga, B. Basturk, Artificial bee colony (abc) optimization algorithm for solving constrained optimization problems, in: P. Melin, O. Castillo, L. T. Aguilar, J. Kacprzyk, W. Pedrycz (Eds.), Foundations of Fuzzy Logic and Soft Computing, Springer Berlin Heidelberg, Berlin, Heidelberg, 2007, pp. 789–798.
  * [59] E.-G. Talbi, Metaheuristics: From Design to Implementation, Vol. 74, John Wiley, 2009.  doi:10.1002/9780470496916.
  * [60] R. W. Erickson, D. Maksimovic, Fundamentals of Power Electronics, 2nd Edition, Springer, 2001.
  * [61] J. Bao, L. Peter, Process Control: The Passive Systems Approach, Springer Publishing Company, Incorporated, 2007.  doi:10.1007/978-1-84628-893-7.
  * [62] M. K. Heris, The Yarpiz Project (2015).
URL https://yarpiz.com/about

  * [63] E. V. Kassick, Deriving the canonical equivalent circuit for small signal & low frequency ac model for the sepic and zeta pwm dc-dc converters with two-port network (quadripoles) circuit analysis technique, Eletr. Potên. 16 (2011) 376–382.  doi:http://dx.doi.org/10.18618/REP.20114.376382.


#  Further stability results of the functional equation
$f(2x+y)+f\left(\frac{x+y}{2}\right)=\frac{2f(x)f(y)}{f(x)+f(y)}+\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}$

I. Sadani  Department of Mathematics, University of Mouloud Mammeri, Tizi-
Ouzou 15000, Algeria<EMAIL_ADDRESS>
###### Abstract.

In this paper, we investigate the generalized Hyers-Ulam stability of the
following reciprocal functional equation

$f(2x+y)+f\left(\frac{x+y}{2}\right)=\frac{2f(x)f(y)}{f(x)+f(y)}+\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}$

in non-Archimedean space using a direct method.

###### Key words and phrases:

Reciprocal Functional Equation, non-Archimedean space, Hyers-Ulam-Rassias
stability.

###### 2010 Mathematics Subject Classification:

39B52, 39B72

##  1\. Introduction

Stability of functional equations have been studied by many authors and is by
now a classical subject. It at least goes back to the question posed by S. M.
Ulam [25] in $1940$. The first affirmative answer to this question was given
by Hyers [5] in $1941$. It is noteworthy that the result of Hyers marked the
inauguration of a series of academic publications in the stability theory of
functional equations. In $1978$, Th. M. Rassias [18] gave a remarkable
generalization of the Hyers’s result which allows the Cauchy difference to be
unbounded. Latterly, The problem of stability for different functional
equations and more in different spaces have been considerably explored by many
authors ( see for example [4],
[1],[3],[2],[7],[6],[8],[9],[10],[11],[12],[13],[16]).

The functional equation

(1.1) $f(x+y)=\frac{f(x)f(y)}{f(x)+f(y)}\ \ \ (x,y\in\mathbb{R}-\\{0\\}).$

is called a reciprocal functional equation. In particular, every solution of
the reciprocal functional equation is said to be a reciprocal mapping. In
$2010$, a generalized Hyers-Ulam stability problem for the reciprocal
functional equation was proved by Ravi and Kumar [14].

Since then, numerous reciprocal functional equations and reciprocal-quadratic,
cubic, quartic, quintic and so on have been discovered and studied by
mathematicians with many interesting stability results (see for example [17,
19, 20, 21, 22, 23, 24] ).

In [19], the author investigated the general solutions and stability of the
following reciprocal type functional equation

(1.2)
$f(2x+y)+f\left(\frac{x+y}{2}\right)=\frac{2f(x)f(y)}{f(x)+f(y)}+\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}$

with $x,y\in X=\mathbb{R}-\left\\{-c\right\\}$ where $c$ is the constant
appears in the general solution

$f(x)=\frac{a}{x+c}.$

In this paper, we investigate generalized Hyers-Ulam-Rassias stability of
functional equation (1.2) in non-Archimedean space.

##  2\. Preliminaries

In this section, we recall some basic definitions and properties about non-
Archimedean fields.

######  Definition 1.

  * •

A valuation is a function $|\cdot|$ from a field $\mathbb{K}$ into
$[0,\infty)$ such that $0$ is the unique element having the $0$ valuation,
$|xy|=|x||y|$ and the following strong triangle inequality holds:

(2.1) $|x+y|\leq\max\\{|x|,|y|\\}.$

  * •

A field $\mathbb{K}$ with a valuation on it is called a valued field.

  * •

A valuation which satisfies the triangle inequality (2.1) is called a non-
Archimedean valuation and the field is called a non-Archimedean field.

######  Remark.

  1. (1)

Let $(\mathbb{K},|\cdot|)$ be a non-Archimedean valued field. Then, for each
$n\in\mathbb{N},|n.1|=|1+1+\cdots+1|\leq\max(|1|,|1|,\cdots,|1|)=1.$

  2. (2)

Let $\mathbb{K}$ be a field. The map defined by $|x|=1$ for $x\neq 0$ and
$|0|=0$ is called trivial valuation and is a non-Archimedean valuation.

######  Example (Non-Archimedean fields).

  1. (1)

For any prime number $p$ the completion $\mathbb{Q}_{p}$ of the field
$\mathbb{Q}$ of rational numbers with respect to the $p-$adic valuation
$|\cdot|_{p}$ given by $|a|_{p}=p^{-r}$ if $a=p^{r}\frac{m}{n}$ and the
integers $m,n$ are not divided by $p$, is a non-archimedean valued field. Its
elements are called $p-$adic numbers.

  2. (2)

Let $p$ be an irreducible polynomial in $\mathbb{R}[x]$. The $p-$adic
valuation $|\cdot|_{p}$ on the rational function field $\mathbb{R}(x)$ is
defined by $|0|_{p}=0$ and $\left|p^{u}\frac{f}{h}\right|=\frac{1}{e^{u}}$
where $u\in\mathbb{Z}$, and $f,h\in\mathbb{R}[x],h\neq 0$ are not divisible by
$p$. This valuation is non-Archimedean and $\mathbb{R}(x)$ is a non-
Archimedean valued field (for more details, see [26]).

######  Definition 2.

Let $Y$ be a vector space over a field $\mathbb{K}$ with a non-Archimedean
valuation $|\cdot|$. A function $\left\|\cdot\right\|:Y\rightarrow[0,\infty)$
is called a non-Archimedean norm if the following conditions hold:

  1. (1)

$\left\|x\right\|=0$ if and only if $x=0$ for all $x\in Y$,

  2. (2)

$\left\|ax\right\|=|a|\left\|x\right\|$ for all $a\in\mathbb{K}$ and $x\in Y$,

  3. (3)

the strong triangle inequality holds:

$\left\|x+y\right\|\leq\max\left\\{\left\|x\right\|,\left\|y\right\|\right\\},$

for all $x,y\in Y$. Then $\left(Y,\left\|\cdot\right\|\right)$ is called a
non-Archimedean normed space.

######  Definition 3.

Let $\\{x_{n}\\}_{n=1}^{\infty}$ be a sequence in a non-Archimedean normed
space $Y.$

  1. (1)

A sequence $\\{x_{n}\\}_{n=1}^{\infty}$ in a non-Archimedean space is a Cauchy
sequence if and only if, the sequence $\\{x_{n+1}-x_{n}\\}_{n=1}^{\infty}$
converges to zero.

  2. (2)

The sequence $\\{x_{n}\\}$ is said to be convergent if, for any $\epsilon>0$,
there are a positive integer $N$ and $x\in Y$ such that

$\left\|x_{n}-x\right\|\leq\epsilon,$

for all $n\geq N$. Then the point $x\in Y$ is called the limit of the sequence
$\\{x_{n}\\}$, which is denote by $\lim_{n\to\infty}x_{n}=x.$

  3. (3)

If every Cauchy sequence in $Y$ converges, then the non-Archimedean normed
space $X$ is called a complete non-Archimedean space.

##  3\. Stability of functional equation (1.2)

Hereinafter, we suppose that $X$ defined as above and $L$ be a complete non-
Archimedean field. We are now in a position to state our main result.

######  Theorem 1.

Let $\mu:X^{2}\rightarrow L$ be a mapping such that

(3.1) $\lim_{n\to\infty}|2|^{n}\mu(2^{n+1}x,2^{n+1}y)=0$

for all $x,y\in X$ and let for each $x\in X$ the limit

(3.2) $\Psi(x)=\lim_{n\rightarrow\infty}\max\\{|2|^{k}\mu(0,2^{k+1}x);0\leq
k<n\\}$

exists. Assume that $f:X\rightarrow L$ is a mapping satisfies

(3.3)
$\left\|f(2x+y)+f\left(\frac{x+y}{2}\right)-\frac{2f(x)f(y)}{f(x)+f(y)}-\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}\right\|\leq\mu(x,y).$

Then

(3.4) $g(x)=\lim_{n\rightarrow\infty}2^{n}f(2^{n}x)$

exists for all $x\in X$ and defines a reciprocal mapping $g:X\rightarrow L$
such that

(3.5) $\left\|f(x)-g(x)\right\|\leq\Psi(x)$

for all $x\in X$. Moreover, if

(3.6) $\lim_{j\rightarrow\infty}\lim_{n\rightarrow\infty}\max\\{|2|^{k}\mu(0,\
2^{k+1}x);j\leq k<n+j\\}=0$

then $g$ is the unique reciprocal mapping satisfying (3.5).

###### Proof.

Substitute $(x,y)$ by $(0,2x)$ in (3.3), we get

(3.7)
$\left\|f\left(x\right)-\frac{2f(0)f(2x)}{f(0)+f(2x)}\right\|\leq\mu(0,2x).$

We replace $x$ by $2^{n}x$ and multiplying by $2^{n}$ in (3.7), we obtain

$\left\|2^{n}f\left(2^{n}x\right)-\frac{2^{n+1}f(0)f(2^{n+1}x)}{f(0)+f(2^{n+1}x)}\right\|\leq|2|^{n}\mu(0,2^{n+1}x)$

which is equivalent to

(3.8)
$\left|2^{n}f\left(2^{n}x\right)-\frac{1}{\frac{1}{2^{n+1}f(2^{n+1}x)}+\frac{1}{2^{n+1}f(0)}}\right|\leq|2|^{n}\mu(0,2^{n+1}x)$

for all $x\in X$. It follows from (3.1), the right-hand side of (3.8) tends to
zero as $n\to\infty$. Thus the sequence $\\{2^{n}f\left(2^{n}x\right)\\}$ is a
Cauchy sequence. Since $L$ is complete, so $\\{2^{n}f\left(2^{n}x\right)\\}$
is convergent. We set

$g(x)=\lim_{n\to\infty}2^{n}f\left(2^{n}x\right).$

By induction, it is straightforward to obtain that

(3.9) $\left\|g(x)-f(x)\right\|\leq\max\\{|2|^{k}\mu(0,2^{k+1}x);0\leq k<n\\}$

for all $n\in\mathbb{N}$ and all $x\in X$. By taking $n$ to approach infinity
in (3.9), and using (3.2), we get (3.5). By (3.1) and (3.3), we have

(3.10)
$\left\|f(2x+y)+f\left(\frac{x+y}{2}\right)-\frac{2f(x)f(y)}{f(x)+f(y)}-\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}\right\|\\\
=\lim_{n\to\infty}||2^{n}f(2^{n+1}(2x+y))+2^{n}f\left(\frac{2^{n+1}(x+y)}{2}\right)-\frac{2^{n+1}f(2^{n+1}x)f(2^{n+1}y)}{f(2^{n+1}x)+f(2^{n+1}y)}\;\\\
-\frac{2^{n+1}f(2^{n+1}(x+y))f(2^{n+1}(y-x))}{3f(2^{n+1}(y-x))-f(2^{n+1}(x+y))}||\leq\lim_{n\to\infty}|2|^{n}\mu(2^{n+1}x,2^{n+1}y)=0.$

for all $x,y\in X$. Therefore the function $f$ : $X\rightarrow L$ satisfies
(1.2). To show the uniqueness of $g$, let $h:X\rightarrow L$ be another
function satisfying (3.5). Then

$\displaystyle\left\|g(x)-h(x)\right\|$
$\displaystyle=\lim_{n\rightarrow\infty}|2|^{n}\left\|g(2^{n+1}x)-h(2^{n+1}x)\right\|$
$\displaystyle\leq\lim_{k\rightarrow\infty}|2|^{n}\max\\{\left\|g(2^{n+1}x)-h(2^{n+1}x)\right\|,\left\|g(2^{n+1}x)-h(2^{n+1}x)\right\|\\}$
$\displaystyle\leq\lim_{j\rightarrow\infty}\lim_{n\rightarrow\infty}\max\\{|2|^{k}\mu(0,2^{k+1}x);j\leq
k<n+j\\}$ $\displaystyle=0$

for all $x\in X$. Therefore $g=h$, and the proof is complete. ∎

The following corollaries are an immediate consequence of Theorem 1.

######  Corollary 1.

Let $\epsilon>0$ be a constant. If $f:X\rightarrow L$ satisfies

(3.11)
$\left\|f(2x+y)+f\left(\frac{x+y}{2}\right)-\frac{2f(x)f(y)}{f(x)+f(y)}-\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}\right\|\leq\epsilon.$

for all $x,y\in X$, then there exists a unique reciprocal mapping $g$ :
$X\rightarrow L$ satisfying (1.2) and $||f\mathcal{}(u)-g(u)||\leq\epsilon$
for all $x\in X.$

###### Proof.

Letting $\mu(x,y)=\epsilon$ in Theorem 1, we obtain the required result. ∎

######  Corollary 2.

Let $\epsilon\geq 0$ and $a\neq-1$, be fixed constants. If $f:X\rightarrow L$
satisfies

(3.12)
$\left\|f(2x+y)+f\left(\frac{x+y}{2}\right)-\frac{2f(x)f(y)}{f(x)+f(y)}-\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}\right\|\\\
\leq\epsilon(|x|^{a}+|y|^{a}).$

for all $x,y\in X$, then there exists a unique reciprocal mapping
$g:X\rightarrow L$ satisfying (1.2) and

$|f(x)-g(x)|\leq\left\\{\begin{array}[]{ll}|2|^{2}\epsilon|x|^{a},&a>-1,\\\
\epsilon\frac{|x|^{a}}{|2|^{a-1}},&a<-1,\end{array}\right.$

for all $x\in X.$

###### Proof.

Considering $\mu(x,\ y)=\epsilon(|x|^{a}+|y|^{a})$ in Theorem 1 the desired
result follows directly. ∎

######  Corollary 3.

Let $\tau:[0,\infty)\rightarrow]0,+\infty]$ be a function satisfying

$\tau(|2|t)\leq\tau(|2|)\tau(t)\ (t\geq 0),\tau(|2|)<|2|.$

Let $\delta>0$ and $f:X\rightarrow L$ be a mapping satisfying

$\left\|f(2x+y)+f\left(\frac{x+y}{2}\right)-\frac{2f(x)f(y)}{f(x)+f(y)}-\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}\right\|\\\
\leq\delta(\tau(|x|)+\tau(|y|))$

with $x,y\in X$. Then there exists a unique reciprocal mapping $h:X\rightarrow
L$ such that

$\|f(x)-h(x)\|\leq\delta\tau(2|x|)\ \ (x\in X).$

###### Proof.

Defining $\mu$ : $X\times X\rightarrow]0,+\infty]$ by $\mu(x,\
y)=\delta(\tau(|x|)+\tau(|y|))$. On the one hand, we have

$\lim_{n\to\infty}|2|^{n}\mu(2^{n+1}x,2^{n+1}y)\leq\lim_{n\to\infty}|2|^{n}\tau(|2|)^{n+1}\mu(x,y)=\lim_{n\to\infty}(|2|\tau(|2|))^{n+1}\mu(x,\
y)$

$=\lim_{n\to\infty}\tau(|2|)(|2|\tau(|2|))^{n}\mu(x,y)=0.$

for all $x,y\in X$. On the other hand,

$\Psi(x)=\lim_{n\rightarrow\infty}\max\\{|2|^{k}\mu(0,2^{k+1}x);0\leq
k<n\\}=\mu(0,2x)=\delta(\tau(|2x|)),$

and

$\lim_{j\rightarrow\infty}\lim_{n\rightarrow\infty}\max\\{|2|^{k}\mu(0,\
2^{k+1}x);j\leq k<n+j\\}=\lim_{j\to\infty}|2|^{j}\mu(0,\ 2^{j+1}x)=0.$

Applying Theorem 1, we get desired result. ∎

##  4\. Conclusion

We have proved the stability of a reciprocal type functional equation

$f(2x+y)+f\left(\frac{x+y}{2}\right)=\frac{2f(x)f(y)}{f(x)+f(y)}+\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}$

in non-Archimedean space using the direct method.

## References

  * [1] B. Belaid, E. Elqorachi, Th. M. Rassias, On the Hyers–Ulam stability of approximately Pexider mappings, Math. Inequal. Appl., 11(2008), 805–818.
  * [2] H. Y. Chu, D. S. Kang, Th. M. Rassias, On the stability of a mixed n-dimensional quadratic functional equation, Bull. Belg. Math. Soc. –Simon Stevin 15(2008) , 9–24.
  * [3] Y. J. Cho, Th. M. Rassias, R. Saadati, Stability of Functional Equations in Random Normed Spaces, Springer, New York, 2013.
  * [4] A. Grabiec, The generalized Hyers-Ulam stability of a class of functional equations, Publ. Math. Debrecen, 48(1996), 217–235.
  * [5] D. H. Hyers. On the stability of the linear functional equation, Proc. Nat. Acad. Sci. U.S.A., 27(1941), 222–224.
  * [6] G. Isac and Th. M. Rassias, Stability of $\psi$-additive mappings: applications to nonlinear analysis, Int. J. Math. Math. Sci., 19(1996), 219–228.
  * [7] G. Isac, Th. M. Rassias, On the Hyers–Ulam stability of additive mappings, J. Approx. Theory, 72(1993), 131–137.
  * [8] S. M. Jung, On the Hyers-Ulam-Rassias st ability of approximately additive mappings, J. Math. Anal. Appl. 204(1996), 221–226.
  * [9] S. M. Jung, Hyers-Ulam-Rassias Stability of Functional Equations in Mathematical Analysis, Hadronic Press, Palm Harbor 2001.
  * [10] S. M. Jung, Hyers–Ulam–Rassias St ability of F unctional Equations in Nonlinear Analysis, Springer Optimization and Its Applications, 48(2011).
  * [11] S. M. Jung, Th. M. Rassias, C. Mortici, On a functional equation of trigonometric type, Appl. Math. Comput. 252(2015), 294–303.
  * [12] M. S. Moslehian, Th. M. Rassias , Orthogonal stability of additive type equations, Aequ. Math. 73(2007), 249–259.
  * [13] A. Najatim, Th. M. Rassias, Stability of a mixed functional equation in several variables on Banach modules. Nonlinear Anal. 72(2010), 1755–1767.
  * [14] K. Ravi, B.V. Senthil Kumar, Ulam-Gavruta-Rassias stability of Rassias reciprocal functional equation. Glob. J. Appl. Math. Math. Sci., 3(2010), 57–79.
  * [15] K. Ravi, J. M. Rassias, B. V. Senthil Kumar, Ulam stability of generalized reciprocal functional equation in several variables, Int. J. Appl. Math. Stat. 19(2010), 1–19.
  * [16] R. Saadati, Y. J. Cho, J. Vahidi, The stability of the quartic functional equation in various spaces, Comput. Math. Appl. 60(2010) , 1994–2002.
  * [17] Ravi, K., Rassias, J.M. and Senthil Kumar, B.V. Ulam stability of a generalized reciprocal type functional equation in non-Archimedean fields. Arab. J. Math. 4, 117–126 (2015).
  * [18] Th.M. Rassias, On the stability of linear mapping in Banach spaces,Proc. Amer. Math. Soc.,72, 297–300 (1978).
  * [19] Sadani, I. On the stability of $f(2x+y)+f\left(\frac{x+y}{2}\right)=\frac{2f(x)f(y)}{f(x)+f(y)}+\frac{2f(x+y)f(y-x)}{3f(y-x)-f(x+y)}$. Tatra Mt. Math. Publ. 78, 1–10 (2020).
  * [20] K. Ravi, J. M. Rassias, B. V. Senthil Kumar, Ulam stability of generalized reciprocal functional equation in several variables, Int. J. Appl. Math. Stat. 19(2010)pp. 119.
  * [21] K. Ravi, J. M. Rassias, B. V. Senthil Kumar, Ulam stability of reciprocal difference and adjoint funtional equations, Australian J. Math. Anal.Appl. 8(1) (2011), pp. 118.
  * [22] K. Ravi, E. Thandapani, B.V. Senthil Kumar, Stability of reciprocal type functional equations, Panam. Math. J. Vol. 21(1) (2011), pp. 59–70.
  * [23] K. Ravi, E. Thandapani, B.V. Senthil Kumar, Solution and stability of a reciprocal type functional equation in several variables, J. Nonlinear Sci. Appl., 7 (2014),pp. 18–27.
  * [24] Kim, S.O., Senthil Kumar, B.V. and Bodaghi, A. Approximation on the reciprocal-cubic and reciprocal-quartic functional equations in non-Archimedean fields, Adv Differ Equ, 77 (2017)
  * [25] S. M. Ulam, Problems in Modern Mathematics, Science Editions, John Wiley and Sons, 1964.
  * [26] Paulo Ribenboim, The theory of classical valuations, Springer Monographs in Mathematics, Springer-Verlag, New York, 1999.