applied
sciences
A icle
A Clus e ing-Based Hyb id Suppo Vec o
Reg ession Model o P edic Con aine Volume a
Seapo Sani a y Facili ies
Juan Jesús Ruiz-Aguila 1,* , JoséAn onio Moscoso-López 1, Daniel U da 2,
Ja ie González-En ique 3and Ignacio Tu ias 3
1Depa men o Indus ial and Ci il Enginee ing, Poly echnic School o Enginee ing, Uni e si y o Cadiz,
11202 Algeci as, Spain; [email p o ec ed]
2G upo de In eligencia Compu acional Aplicada (GICAP), Depa amen o de Ingenie ía In o má ica,
Escuela Poli
é
cnica Supe io , Uni e sidad de Bu gos, A . Can ab ia s/n, 09006 Bu gos, Spain; du [email p o ec ed]
3Depa men o Compu e Science Enginee ing, Poly echnic School o Enginee ing, Uni e si y o Cadiz,
11202 Algeci as, Spain; ja ie [email p o ec ed] (J.G.-E.); [email p o ec ed] (I.T.)
*Co espondence: juanjesus. [email p o ec ed]
Recei ed: 7 No embe 2020; Accep ed: 22 No embe 2020; Published: 24 No embe 2020


Abs ac :
An accu a e p edic ion o eigh olume a he sani a y acili ies o seapo s is a key ac o
o imp o e planning ope a ions and esou ce alloca ion. This s udy p oposes a hyb id app oach
o o ecas con aine olume a he sani a y acili ies o a seapo . The me hodology consis s o a
h ee-s ep p ocedu e, combining he s eng hs o linea and non-linea models and he capabili y
o a clus e ing echnique. Fi s , a sel -o ganizing map (SOM) is used o decompose he ime se ies
in o smalle clus e s easie o p edic . Second, a seasonal au o eg essi e in eg a ed mo ing a e ages
(SARIMA) model is applied in each clus e in o de o ob ain p edic ed alues and esiduals o
each clus e . These alues a e inally used as inpu s o a suppo ec o eg ession (SVR) model
oge he wi h he his o ical da a o he clus e . The inal p edic ion esul in eg a es he p edic ion
esul s o each clus e . The expe imen al esul s showed ha he p oposed model p o ided accu a e
p edic ion esul s and ou pe o ms he es o he models es ed. The p oposed model can be used as
an au oma ic decision-making ool by seapo managemen due o i s capaci y o plan esou ces in
ad ance, a oiding conges ion and ime delays.
Keywo ds:
ma i ime anspo ; con aine o ecas ing; suppo ec o eg ession; sel -o ganizing maps;
machine lea ning; hyb id models
1. In oduc ion
O e he las decades, po s ha e played an impo an ole in in e na ional ade and mos o he
o e seas shipping o p oduc s is aboa d deep-sea con aine essels [
1
]. The inc ease in a ic o goods
and he Eu opean uni ica ion has led o conside he enhancemen o secu i y a bo de c ossings
o he Eu opean Union. In his sense, he Bo de Inspec ion Pos s (BIPs) we e c ea ed in o de o
gua an ee he secu i y a bo de c ossings and he quali y o he impo -expo goods. BIPs a e he
app o ed acili ies whe e he checks o goods ( anspo ed wi hin con aine s by ucks o owing
ehicles) a e ca ied ou be o e en e ing he Communi y e i o y. Howe e , he sus ained g ow h in
he wo ldwide exchange o goods is c ea ing he need o u he inspec ions esul ing in conges ion
and high load-peaks wi hin he sani a y acili ies. This causes ime delays and highe cos s in he
supply chain. The BIPs a e he eby bo lenecks ha mus be necessa ily assessed by Po Au ho i ies in
o de o keep he quali y le el o he po and a oid losing compe i i eness. In o de o a oid ime
delays and conges ion in he sani a y acili ies, he po managemen should be able o accu a ely
Appl. Sci. 2020,10, 8326; doi:10.3390/app10238326 www.mdpi.com/jou nal/applsci
Appl. Sci. 2020,10, 8326 2 o 17
o ecas he numbe o con aine passing h ough hese acili ies. An accu a e p edic ion o he olume
o con aine a ic h ough he BIP may become a use ul ool o imp o e human esou ces, planning
ope a ions and he se ice quali y a po s.
Fo ecas ing ime se ies ha e ocused a en ion o many esea che s o a long ime. Many e o s
ha e been ca ied ou o imp o e he exis ing me hodologies and o achie e enhanced models o
ob ain accu a e p edic ions in any o ecas ing ields ( o ins ance, ene gy consump ion, anspo a ion,
en i onmen , economy, medicine and heal h). The ocus o his s udy is mainly se on ma i ime
anspo due o he a ic associa ed wi h he po BIPs. The p oposed o ecas ing echniques can be
di ided in o h ee ca ego ies: single me hods, combined me hods and hyb id me hods.
The i s class comp ises bo h linea and nonlinea echniques. Linea echniques a e based on he
assump ion ha a linea ela ionship exis s be ween he u u e alues and he cu en and pas alues
o he ime se ies. One o he linea models ha has a ac ed mo e a en ion o e he pas ew decades
is he well-known au o eg essi e in eg a ed mo ing a e ages (ARIMA) model and i s u he ex ended
e sions. Based on he Box and Jenkins me hodology [
2
] he ARIMA models ha e been cons an ly
applied o sol e o ecas ing asks ela ed o ma i ime anspo . ARIMA models we e success ully
applied o p edic expo ope a ions in con aine e minals [
3
], o p edic ce ain a ic lows o goods
ha pass h ough a po [4,5] and ecen ly o p edic ing con aine h oughpu olumes a po s [6].
Mo eo e , nonlinea echniques ha e become a s eng h al e na i e agains he weaknesses o
linea models. This o ecas ing echnique has p o ed o be e ec i e when ime se ies show nonlinea
pa e ns, o e coming he main cons ain s o linea models. Real-wo ld p oblems a e complex and
mos ly gene a ed by an unde lying nonlinea p ocess [
7
], such could be he case o olume o con aine s
demand a BIPs. In his subca ego y, wo machine lea ning echniques highligh : a i icial neu al
ne wo ks (ANNs) and suppo ec o machines o eg ession (SVR). Ins ead o ANNs, SVR has
a ac ed inc easing a en ion in ecen yea s due o i s inhe en abili ies o o e come some gene al
limi a ions o ANNs, such as he achie emen o a global minimum. Due o i s g ea gene aliza ion
abili y, SVR has been used in o ecas ing anspo asks wi h p omising esul s. These wo echniques
ha e been cons an ly compa ed in he esea ch li e a u e.
The second ca ego y comp ises he combined models. These models a emp o o e come he
cons ain s o he single o ecas ing models when ime se ies exhibi seasonal ea u es, empo al
e olu ions, co ela ions, and o he pa e ns ha d o cap u e by single models. To his aim, single
me hods a e combined o seize he abili ies o each one in ce ain o ecas ing asks. One o he
mos equen ly used app oaches consis s o combining a single p edic ion echnique (mainly a so
compu ing echnique) wi h a clus e ing me hod. Based on he di ide-and-conque p inciple, clus e ing
me hods allow o di ide he o iginal da abase in o a numbe o smalle g oups, called clus e s. The main
assump ion is ha simple clus e s a e easie o sol e han add essing he whole da abase. When he
clus e ing me hod has di ided he da abase in o se e al clus e s, a p edic ion echnique is hen applied
in each clus e independen ly. Sel -o ganizing maps (SOMs) [
8
] is p obably he bes -known clus e ing
me hod. This popula unsupe ised lea ning echnique has been p o ed o be as e and mo e accu a e
and e icien han o he clus e ing me hods [
9
], imp o ing he inal p edic ion pe o mance in ime
se ies [
10
]. A combined SOM-ANN model was in oduced by Chen e al. [
11
] o p edic a ic lows in
anspo a ion. This wo k also compa ed he o ecas ing pe o mance o he p oposed model o he
ob ained wi h a SOM-ARIMA model and a single ARIMA model. Resul s showed ha he SOM-ANN
model ou pe o med he es o models. Due o he ecen eme gence o SVR in anspo a ion,
he e is ha dly any esea ch ela ed o anspo combining SOM and SVR in a wo-s age p ocedu e.
Ne e heless, i is a widesp ead solu ion in many o he o ecas ing ields [12].
The hi d ca ego y includes hyb id models. Single linea models ha e shown o ha e se e al
cons ain s, limi a ions and disad an ages in pa icula si ua ions and ce ain applica ions ( o ins ance,
when he unde lying gene a ing mechanism is nonlinea o unce ain y is p esen ed) [
13
]. In a simila
manne , he use o nonlinea models can be o ally inapp op ia e i linea pa e ns a e p esen ed.
Mo eo e , eal-wo ld ime se ies a e no comple ely linea o nonlinea , bu a he con ain bo h
Appl. Sci. 2020,10, 8326 3 o 17
componen s. Thus, a me hodology using linea and nonlinea models in a hyb id way akes he
capabili ies o bo h models, imp o ing he o ecas ing accu acy and educing he isk o ailu e when
an unsui able single model is used. Hyb idizing linea models and machine lea ning echniques ha e
been p oposed in ecen yea s o o ecas anspo a ion ime se ies. Pa icula ly, ARIMA has been he
mos commonly used linea model in li e a u e o cons uc his kind o hyb id model. The i s wo ks
we e p oposed conside ing ANNs as he machine lea ning echnique in [
7
,
14
] o o ecas ime se ies
and hey concluded ha hyb id ARIMA-ANN models we e supe io o single SARIMA and ANN
models. Since hen, many esea ch wo ks can be ound in he li e a u e in many a eas linked o ime
se ies analysis [
15
,
16
]. Se e al au ho s we e also p oposed a hyb idiza ion o SARIMA and SVR o
add ess se e al o ecas ing asks ou side he anspo sec o [
17
,
18
], al hough i is less widesp ead
compa ed o SARIMA-ANN models. Rela ed o ma i ime anspo , Xie e al. [
19
] p oposed se e al
hyb id app oaches in a compa a i e way including he SARIMA-SVR model o con aine h oughpu
o ecas ing. All hese p e ious s udies coincided in poin ing ou ha a hyb id s a egy conside ing
ARIMA and SVR models o e came he pe o mance o single models in hei espec i e domains.
In his s udy, a combined-hyb id o ecas ing model is p oposed in such a way ha a hyb id model
(SARIMA-SVR) is combined wi h a clus e ing me hod (SOM) o o ecas he daily numbe o con aine s
passing h ough a BIP, hus esul ing in a SOM-SARIMA-SVR model. This me hodology uni ies he
s eng hs o clus e ing me hods in decomposing he o ecas ing ask in o some ela i ely easie sub asks
(using a SOM me hod) and he s eng hs o hyb id models o i linea and nonlinea componen s
(using a SARIMA-SVR model). Thus, he inal aim o his s udy is wo old: i s , o demons a e ha he
SOM-SARIMA-SVR ou pe o ms he es o possible hyb id o combined models in o ecas ing he daily
numbe o con aine s passing h ough BIP o a ma i ime po ; and second, o es he gene aliza ion
capabili ies o SVR in o ecas ing logis ic asks, especially he con aine inspec ion p ocess ha can be a
bo leneck in he supply chain. To ain he new model, a h ee-s ep p ocedu e was de eloped. In he
i s s ep a SOM algo i hm allows o di ide he da abase in o di e en clus e s o egions. On a second
s ep, a SARIMA model is i ed o he da a o each clus e in o de o cap u e he seasonali y and he
linea beha io . Finally, a SVR model is applied o e each ob ained clus e conside ing se e al hyb id
app oaches (di e en inpu con igu a ions) which can conside he o iginal da abase and he ou pu s
o he i s s ep ( esidual and p edic ed alues).
The es o he pape is o ganized as ollows: The second sec ion gi es an o e iew o
sel -o ganizing maps (SOM), seasonal au o eg essi e in eg a ed mo ing a e age (SARIMA) and
suppo ec o machines o eg ession (SVR). The empi ical da a, he combined-hyb id me hodology,
and he expe imen al p ocedu e a e p esen ed in Sec ion 3. The ou h sec ion discusses he esul s
ob ained. Finally, he las sec ion summa izes he impo an conclusions o his wo k.
2. Me hods
The p oposed model consis s o hyb idizing a seasonal au o eg essi e in eg a ing mo ing a e age
(SARIMA) model as a linea model, and a non-linea model such as a suppo ec o machine
o eg ession (SVR). These me hods a e wo o he mos impo an o ecas ing models conce ning
hei espec i e domains. Addi ionally, a clus e ing me hod, Kohonen sel -o ganizing map (SOM),
is combined wi h his hyb id (SARIMA-SVR) model. The me hodology and basic concep s a e
nex desc ibed.
2.1. Sel -O ganizing Maps (SOM)
Wi hin he unsupe ised lea ning ield, a SOM is a kind o neu al ne wo k ha has gained special
a en ion du ing he pas wo decades. Fi s p oposed by Kohonen [
8
,
20
], a SOM is a classi ica ion
echnique ha g oups objec s o he sys ems in o egions called clus e s. This classi ica ion is based
on he simila i y o nea ness o hese objec s wi hou ex e nal ac o s in luencing hei pe o mance.
In he p ocess, he neu ons o he model o ganize hemsel es conside ing only hose ha play a simila
Appl. Sci. 2020,10, 8326 4 o 17
ole, o ming a clus e . The ep esen a i e poin o his clus e is called cen oid, being he cen al
poin , which can be used as he cen e poin o a classi ie based on minimum dis ance.
The opology o a SOM model consis s o se e al neu ons dis ibu ed in o wo laye s: he inpu
and he ou pu laye . The i s one is o med by kneu ons. Each neu on co esponds wi h one
inpu . The ou pu laye , called he compe i ion laye , can consis o di e en opologies (2-D g id
o his case). The p e-p ocessing is pe o med in his laye . In he p ocess, all he neu ons jo
he ou pu laye a e connec ed by weigh s (w
i,j
) wi h all neu ons io he inpu laye . A weigh
ec o ,
wj=wj,1,wj,2,. . . ,wj,k
, is he eby associa ed o each ou pu neu on jo he compe i i e laye
(kis he eby he inpu ec o dimension). The aining o a SOM model is based on a compe i i e p ocess
whe e ne wo ks a e ained i e a i ely. Di e en inpu ec o s
xi= [x1,x2,. . . ,xk]T
a e p esen ed o he
ne wo k a each i e a ion. Du ing he ne wo k aining, he Euclidean dis ance be ween xand all he
weigh ec o s a e compu ed as ollows:


x−wb

=min
jnkx( )−wj( )ko,j=1, 2, . . . ,l, (1)
whe e lis he numbe o ou pu neu ons. Acco ding o Equa ion (1), w
b
is conside ed he winning neu on,
i.e., he neu on ha has he weigh ec o closes o x. In addi ion, he weigh o he winning neu on is
upda ed a ime +1, acco ding o Equa ion (2). The same occu s wi h i s associa ed neighbo neu ons:
wi( +1) = βi·δb,i( )·[x( )−wi( )]+wi( )(2)
whe e
δb,i
( ) is he neighbo hood unc ion (usually a Gaussian unc ion) associa ed o he neu on i,
and
β
( ) is he exponen ial decay lea ning ac o . Bo h pa ame e s,
β
( ) and
δb,i
( ), decay wi h ime.
The aining algo i hm s ops when he maximum numbe o epochs (a ep esen a ion o all inpu s
pa e ns) is achie ed, o when he pe o mance is minimized o he a ge .
2.2. Au o-Reg essi e In eg a ed Mo ing A e ages (ARIMA)
ARIMA models we e in oduced by Box and Jenkins [
2
]. ARIMA has been a widely used
o ecas ing linea model du ing se e al decades. In hese kind o models, he u u e alue is accep ed o
be a linea unc ion o se e al pas obse a ion and an e o e m. Th ee p edic ion e ms compose his
linea unc ion: he au o eg essi e e m (AR), he mo ing a e age e m (MA), and he in eg a ion e m
(I). A SARIMA model can be ob ained by ex ending he ARIMA model o include seasonal ea u es.
In his way, he model is speci ied as SARIMA(p,d,q)(P,D,Q)
S
, whe e q ep esen he o de o he mo ing
a e age e ms, pdeno es he o de o he au o eg essi e e ms, and dis he deg ee o di e encing.
(P,D,Q) deals wi h he seasonal pa and he capi al le e s co esponds o hei coun e pa s o he
seasonal models wi h he seasonal o de s and he seasonali y o he model is ep esen ed by he
pa ame e s. Equa ion (3) depic s a ypical exp ession o he SARIMA model:
ϕp(L)ΦP(Bs)∇d∇D
sy =θq(B)ΘQ(Bs)a , (3)
whe e y
is he obse ed alue,
∇d
and
∇D
s
a e he egula and seasonal di e encing ope a o s,
espec i ely, pand Pa e he numbe o non-seasonal and seasonal au o eg essi e e ms, qand Q
a e he numbe o non-seasonal and seasonal mo ing a e age e ms, dand Da e he numbe o
egula and seasonal di e ences,
ϕ
and
Φ
depic he alue weigh s o he non-seasonal and seasonal
au o eg essi e e m,
θ
and
Θ
ep esen he weigh s o he non-seasonal and seasonal mo ing a e age
e m, he seasonali y is ep esen ed by S, and a
is he noise e m. The i s s ep is o iden i y he SARIMA
s uc u e assis ed by he obse a ion o he simple and pa ial au oco ela ion unc ion (ACF and
PACF) o he ime se ies. In addi ion, da a mus be s a iona i y in mean, using powe ans o ma ions,
and s a iona i y in a iance, employing di e encing o he ime se ies. Second, he pa ame e s o he
model a e es ima ed. Thi d, he es ima ed esiduals mus be checked. The equi emen s o a whi e
noise p ocess should be sa is ied by he esiduals and hey a e assumed o be independen . Se e al
Appl. Sci. 2020,10, 8326 5 o 17
s a is ic es s and plo s o he esiduals a e used o his pu poses. Finally, an es ima ion o he u u e
alue o he olume o con aine s is ob ained. In his wo k, a alue ange o he model pa ame e s a e
i e a i ely es ed o iden i y he mos sui able model.
2.3. Suppo Vec o Machines o Reg ession (SVR) Models
In con as o ANN models, suppo ec o machines (SVM) is a kind o machine lea ning
echnique ocused on he s uc u al isk minimiza ion ins ead o he empi ical isk minimiza ion
p inciple. The main objec i e o his me hod is maximizing he ma gin dis ance [
21
]. The model can be
o mula ed as he ollowing equa ion:
y(x) = wTφ(x) + b, (4)
whe e bdeno es he bias e m and wis he ec o o weigh s.
φ
(x) ep esen s he ke nel unc ion used
o deal wi h he nonlinea p oblem, mapping he inpu da a in o a highe dimensional ( ea u e) space
whe e da a can be linea . Fi s in oduced o classi ica ion p oblems, he
ε
-insensi i e loss unc ion,
p esen ed in Equa ion (5), has enabled i s use in eg ession p oblems:
Lε(y)(0, i  (x)−y≤ε
 (x)−y−ε, o he wise , (5)
In Equa ion (5),
ε
deno es he a ea o
ε
-insensi i e. The p ocess is he ollowing: i s , he inpu
da a a e mapped in o a new space o highe dimensional ea u es, called ea u e space, by a non-linea
mapping a p io i using a ke nel ans o ma ion. The aim o his ea u e space is o de ec a linea
eg ession unc ion ha can be i he ou pu da a wi h he inpu da a. This linea eg ession co esponds
o he nonlinea eg ession model in he o iginal space. Equa ion (6) ep esen s he p oblem ha should
be op imized:
min
w,b,ξ
1
2||w||2+C
N
P
i=1
(ξ+
i+ξ−
i)
subjec o :
w·xi+b−yi≤ε+ξ+
i
yi−w·xi−bi≤ε+ξ−
i
ξ+
i,ξ−
i≥0
(6)
wi h
i=
1,
. . .
,
l
.
ξi−
and
ξi+
a e he slack a iables ha deal wi h he aining e o on he op
and he bo om, espec i ely. These a iables ake ze o alues wi hin he {
−ε
,
ε
} a ea, and non-ze o
alues ou side i .
1
2||w||2
is he s uc u e isk conce ning he la ness o he model and he pa ame e
Cis a co ec ion ac o ha deals wi h he ade-o be ween he la ness and he e o . In his wo k,
he Gaussian ke nel was chosen as ke nel unc ion. The dual op imiza ion p oblem can be sol ed wi h
he Lag angian mul iplie me hod [
21
]. Finally, he me hod ob ains he suppo ec o s as he inal
decision, which a e he obse a ions wi h non-ze o coe icien s o he Lag angian mul iplie s.
3. Fo ecas ing App oach
In his s udy, a h ee-s ep p ocedu e based on a SOM-SARIMA-SVR model is p oposed in o de
o p edic he daily numbe o con aine s passing h ough a sani a y acili y o a seapo . A u he
objec i e is o compa e he p edic ion esul s o he p oposed h ee-s ep p ocedu e wi h he ob ained
using o he possible me hodologies, such as he SVR models in a single way, he hyb id SARIMA-SVR
model and he combined SOM-SVR model. The p oposed app oach combines he capabili y o SARIMA
models in cap u ing he linea pa e ns and he abili y o SVR in modelling nonlinea pa e ns, oge he
wi h he ad an ages o he SOM clus e ing algo i hm.
The expe imen al da abase comes om he BIP o he Po o Algeci as Bay, loca ed in he Sou h
o Spain. The BIP o Algeci as is he app o ed acili y whe e impo goods a e inspec ed and checked

Appl. Sci. 2020,10, 8326 6 o 17
be o e en e ing he Communi y e i o y o he Eu opean Union. The Po o Algeci as Bay was he
i s po in he Medi e anean Sea and he ou h po in he Eu opean con inen ela ed o he o al
h oughpu in 2019 (109.4 million ons) only su passed by he po s o Ro e dam (The Ne he lands),
An we p (Belgium) and Hambu g (Ge many). The da abase was p o ided by he Po Au ho i y and
con ains daily eco ds o he numbe o con aine s a he Algeci as BIP om 2010 o 2014.
3.1. The P oposed Hyb id Me hodology
A SOM algo i hm (s ep I) was i s ly used o di ide he whole da abase in o se e al disjoin
clus e s wi h simila s a is ical p ope ies. Then, he seasonali y and he linea pa e ns we e cap u ed
using a SARIMA model in each clus e (s ep II). As a esul , a se o p edic ed alues and esiduals
om each clus e was ob ained. Finally, a SVR model was hen implemen ed in he hi d s ep o e
each clus e o gene alize he non-linea ela ionship, ob aining he inal p edic ion (s ep III). Inpu s o
he SVR model we e he o iginal and p edic ed da a om he second (SARIMA) s ep. Two hyb id
app oaches we e es ed in he las s ep conside ing he con igu a ion and he numbe o inpu s.
The o ecas ing pe o mance was assessed using di e en p edic ion ho izons: o he SARIMA
s ep, only one-day (ph =1) p edic ion ho izon was conside ed; and o he inal (SVR) s ep, wo p edic ion
ho izons we e es ed, one-day (ph =1) and se en-day (ph =7) ahead. The p edic ion is one-s ep ahead
(y
+ph
) in bo h cases. In he ph =1 case, he p edic ion is y
+1
, and o he ph =7 case he p edic ion is
y
+7
. The es ima ion can be he eby modelled as a nonlinea unc ion o he np eceding alues o he
ime se ies and an e o e m. This is called he au o eg essi e window (n) and i s design is p esen ed
in Figu e 1. The au o eg essi e window o ph =1 (s eps II and III) is p esen ed on he op o he igu e
and o ph =7 (s ep III) is showed on he bo om o he igu e.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 o 17
The expe imen al da abase comes om he BIP o he Po o Algeci as Bay, loca ed in he Sou h
o Spain. The BIP o Algeci as is he app o ed acili y whe e impo goods a e inspec ed and checked
be o e en e ing he Communi y e i o y o he Eu opean Union. The Po o Algeci as Bay was he
i s po in he Medi e anean Sea and he ou h po in he Eu opean con inen ela ed o he o al
h oughpu in 2019 (109.4 million ons) only su passed by he po s o Ro e dam (The Ne he lands),
An we p (Belgium) and Hambu g (Ge many). The da abase was p o ided by he Po Au ho i y and
con ains daily eco ds o he numbe o con aine s a he Algeci as BIP om 2010 o 2014.
3.1. The P oposed Hyb id Me hodology
A SOM algo i hm (s ep I) was i s ly used o di ide he whole da abase in o se e al disjoin
clus e s wi h simila s a is ical p ope ies. Then, he seasonali y and he linea pa e ns we e cap u ed
using a SARIMA model in each clus e (s ep II). As a esul , a se o p edic ed alues and esiduals
om each clus e was ob ained. Finally, a SVR model was hen implemen ed in he hi d s ep o e
each clus e o gene alize he non-linea ela ionship, ob aining he inal p edic ion (s ep III). Inpu s
o he SVR model we e he o iginal and p edic ed da a om he second (SARIMA) s ep. Two hyb id
app oaches we e es ed in he las s ep conside ing he con igu a ion and he numbe o inpu s.
The o ecas ing pe o mance was assessed using di e en p edic ion ho izons: o he SARIMA
s ep, only one-day (ph = 1) p edic ion ho izon was conside ed; and o he inal (SVR) s ep, wo
p edic ion ho izons we e es ed, one-day (ph = 1) and se en-day (ph = 7) ahead. The p edic ion is one-
s ep ahead (y +ph) in bo h cases. In he ph = 1 case, he p edic ion is y +1, and o he ph = 7 case he
p edic ion is y +7. The es ima ion can be he eby modelled as a nonlinea unc ion o he n p eceding
alues o he ime se ies and an e o e m. This is called he au o eg essi e window (n) and i s design
is p esen ed in Figu e 1. The au o eg essi e window o ph = 1 (s eps II and III) is p esen ed on he
op o he igu e and o ph = 7 (s ep III) is showed on he bo om o he igu e.
Figu e 1. Au o eg essi e window sizes in S eps II and III and hei p edic ion ho izons (ph): one-day
p edic ion ho izon (below he imeline) and se en-day p edic ion ho izon (allow he imeline). n is
he size o he au o eg essi e window.
Fo he ph = 7 case, he au o eg essi e window is composed o alues o he con aine se ies
pe iodically sampled e e y se en days in he pas . This is due o he weekly seasonali y ound in he
analysis o he au oco ela ion unc ion o he ime se ies. The main assump ion he e is ha bes
p edic ion is ob ained when using as inpu s se e al samples om he pas o he same day o he
week (i.e., using se e al successi e Mondays in he pas o p edic a u u e Monday). The p oposed
SOM-SARIMA-SVR p ocedu e is shown g aphically in Figu e 2.
ph = 7
…
n = 1
nn = 2
……
…
y
y -7
y -14
y -n·ph
y
y +7
ph = 1
…
n
y -1
y -2
y -3
y -4
y -5
y -6
y -7
y -n·ph
n = 1n = 2
n = 3
y +1 y +7
…
Pas Fu u e
P esen
Figu e 1.
Au o eg essi e window sizes in S eps II and III and hei p edic ion ho izons (ph): one-day
p edic ion ho izon (below he imeline) and se en-day p edic ion ho izon (allow he imeline). nis he
size o he au o eg essi e window.
Fo he ph =7 case, he au o eg essi e window is composed o alues o he con aine se ies
pe iodically sampled e e y se en days in he pas . This is due o he weekly seasonali y ound in
he analysis o he au oco ela ion unc ion o he ime se ies. The main assump ion he e is ha bes
p edic ion is ob ained when using as inpu s se e al samples om he pas o he same day o he
week (i.e., using se e al successi e Mondays in he pas o p edic a u u e Monday). The p oposed
SOM-SARIMA-SVR p ocedu e is shown g aphically in Figu e 2.
Appl. Sci. 2020,10, 8326 7 o 17
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 o 17
Con aine olume da a
Y
Clus e 2Clus e 1 Clus e n
SOM
y1y2yn
SARIMA 1 SARIMA 2 SARIMA n
e1, p1e2, p2en, pn
SVR 2 SVR n
Hyb id
app oach 1
Hyb id
app oach 2
1 1 1 1
(, , )y pye
11() ee
1 1 1
ype
n
y
2
y
1
y
Final Ou comes
12
, ,...,


n
Yy y y
SVR 1
Hyb id
app oach 1
Hyb id
app oach 2
1 1 1 1
(, , )y pye
11() ee
1 1 1
ype
Hyb id
app oach 1
Hyb id
app oach 2
1 1 1 1
(, , )y pye
11() ee
1 1 1
ype
Figu e 2. The o e all p ocess scheme o he SOM-SARIMA-SVR app oach.
3.1.1. S ep I: SOM
A sel -o ganizing ea u e map (SOM) is i s applied o he da a in o de o spli he da abase
in o se e al disjoin g oups, called clus e s, wi h simila s a is ical dis ibu ion. The assump ion is
ha , in he nex s ep, a o ecas ing echnique can p edic mo e accu a e each clus e ins ead o he
whole da abase. Each clus e wo ks independen ly in he second and hi d s ep. In such cases, a single
SARIMA and SVR models a e applied independen ly a e decomposing he he e ogeneous da a in o
smalle homogeneous egions. A p ede ined numbe o clus e s, c, a e selec ed, be o e s a ing he
SOM algo i hm, o assess and compa e he pe o mance o di e en solu ions.
An expe imen al amewo k was de eloped in o de o selec he op imal numbe o pas alues
(nc) o he inpu ec o ( he numbe o inpu neu ons). Thus, a da abase o inpu s was designed as in
Equa ion (7):
 
12
, , ,..., , T
i k k nc ph
x y y y y
     

,
(7)
whe e is each sample in he daily ime se ies om 2010 o 2014 and k means he sample pe iod
explained g aphically in Figu e 1 (k = 1 o k = 7 days). Each ow is a anged ecu si ely using di e en
lagged e ms (as in an au o eg essi e window).
A la ge numbe o di e en expe imen s we e pe o med modi ying he inpu ec o size and
he numbe o i e a ions. In o de o conside he inhe en andomness o he p ocess, 20 epe i ions
o each expe imen we e ca ied ou .
Figu e 2. The o e all p ocess scheme o he SOM-SARIMA-SVR app oach.
3.1.1. S ep I: SOM
A sel -o ganizing ea u e map (SOM) is i s applied o he da a in o de o spli he da abase
in o se e al disjoin g oups, called clus e s, wi h simila s a is ical dis ibu ion. The assump ion is
ha , in he nex s ep, a o ecas ing echnique can p edic mo e accu a e each clus e ins ead o he
whole da abase. Each clus e wo ks independen ly in he second and hi d s ep. In such cases, a single
SARIMA and SVR models a e applied independen ly a e decomposing he he e ogeneous da a in o
smalle homogeneous egions. A p ede ined numbe o clus e s, c, a e selec ed, be o e s a ing he
SOM algo i hm, o assess and compa e he pe o mance o di e en solu ions.
An expe imen al amewo k was de eloped in o de o selec he op imal numbe o pas alues
(nc) o he inpu ec o ( he numbe o inpu neu ons). Thus, a da abase o inpu s was designed as in
Equa ion (7):
xi=hy ,y −1·k,y −2·k,. . . ,y −nc·ph,iT, (7)
whe e is each sample in he daily ime se ies om 2010 o 2014 and kmeans he sample pe iod
explained g aphically in Figu e 1(k=1 o k=7 days). Each ow is a anged ecu si ely using di e en
lagged e ms (as in an au o eg essi e window).
A la ge numbe o di e en expe imen s we e pe o med modi ying he inpu ec o size and he
numbe o i e a ions. In o de o conside he inhe en andomness o he p ocess, 20 epe i ions o
each expe imen we e ca ied ou .
Appl. Sci. 2020,10, 8326 8 o 17
3.1.2. S ep II: SARIMA
A SARIMA model is i ed o each clus e gene a ed by he SOM in S ep I, ob aining di e en
p edic ed and esidual alues o hese clus e s. The main eason o selec SARIMA as a o ecas ing
model is i s abili y o cap u e linea pa e ns. In o de o assess he o ecas ing pe o mance, and due
o he de e minis ic na u e o linea models, a hold-ou alida ion echnique was applied du ing he
p ocess. The da a (o each clus e ) was di ided in o wo g oups: he aining se con aining wo- hi ds
o he da ase , and he es se comp ising he es o he samples. The pa ame e s o he model we e
adjus ed using he aining se and he es se was used o alida e he model. Di e en alues o he
model pa ame e s we e es ed using a ial-and-e o p ocedu e. Only he ph =1 p edic ion ho izon
was es ed. The alues o he pa ame e s es ed wi hin each clus e a e, o he non-seasonal pa :
p=0, 1, 2, 3, 4;
d=0, 1, 2 and q=1, 2, 3, 4; and o he seasonal pa : s=2, 5, 7; P=0, 1, 2; D=0, 1, 2
and Q=0, 1, 2, 3.
3.1.3. S ep III: SVR
A his s age, a SVR model is applied o each clus e . The h ee di e en adjus able pa ame e s,
which go e ns he SVR model, we e de e mined by an i e a i e p ocess ( ial-and-e o ). The pa ame e s
o he di e en SVR models es ed a e shown in Table 1. No e ha he op imum numbe o clus e s
was wo. This issue is u he explained below.
Table 1.
Pa ame e anges used o he SVR models in he hi d s ep. The subsc ip indica es he clus e
o which he pa ame e belongs o (1 o 2).
Pa ame e s Clus e Values
Hyb id con igu a ion: SOM-SARIMA-SVR-1 SOM-SARIMA- SVR-2
Numbe o a iables:
Inpu laye 1 e1y1,p1,e1
2e2y2,p2,e2
Ou pu laye 1 e1y1
2e2y2
Au o eg essi e window size (s ep: 1):
Residual (e): ne 1, 2 1:20 none, 1:20
O iginal da a (y): ny 1, 2 - none, 1:10
P edic ed alues (p): np 1, 2 - none,1:20
Pa ame e s:
ε(s ep: 1): 1, 2 2(−12:−2) 2(−12:−2)
γ(s ep: 1): 1, 2 2(−12:−2) 2(−12:−2)
C: 1, 2 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100:1000 (s ep: 100)
P edic ion ho izons 1, 2 1, 7 days
Fo each clus e , he inpu s o he SVR model a e o med by he clus e ed da a o he o iginal
ime se ies, y
i
, and hei o ecas ed alues, p
i
, and esiduals, e
i
, ob ained in he second (SARIMA)
s ep (ideno es he clus e hey belong o). The p esence o absence o hese a iables in he inpu s
de e mine he hyb id con igu a ions. Each a iable is so ed ecu si ely as an au o eg essi e window.
The sizes o he o iginal da a, p edic ed alues and esiduals om he SARIMA s ep a e deno ed as ny,
np and ne, espec i ely.
Due o he educ ion o a ailable da a in he clus e s, a andomized esampling s a egy was
implemen ed du ing he SVR p ocess o ensu e he independence o he esul s. Fo his wo k, a wo old
c oss- alida ion (2-CV) echnique was used. This alida ion s a egy was epea ed 20 imes and he
inal alue o he p edic ion pe o mances was he a e age o hese epe i ions. The 20 epe i ions we e
implemen ed o each combina ion o pa ame e alues o he ange shown in Table 1. Once all he ange
alues o he pa ame e s we e es ed, he pa ame e combina ion ha achie es he bes pe o mance
Appl. Sci. 2020,10, 8326 9 o 17
index alues is selec ed. The o al p edic ed ime se ies is o med by adding he p edic ions o each
clus e s. No e ha , as in he SARIMA model, he bes SVR model may be di e en on each clus e .
Two hyb id app oaches we e p oposed and assessed. The p edic ion esul s we e ob ained o
wo p edic ion ho izons, ph =1 and ph =7.
SOM-SARIMA-SVR-1 model (hyb id app oach 1). This classical app oach conside s ha he
ela ionship be ween he linea componen L
and he nonlinea componen NL
o a ime se ies is
addi i e. The e o e, he ime se ies can be decomposed in hese wo independen e ms. Then, a linea
o ecas ing model such as SARIMA can be applied in o de o model he linea componen and he eby
o ob ain he p edic ed alues,
ˆ
p
, and he esidual, e
. Subsequen ly, a SVR model is applied o e
he esiduals in o de o i he nonlinea componen NL
. The main assump ion o his adi ional
app oach is ha he nonlinea ela ionship o he ime se ies can be ound, i any, in he esiduals o he
linea model. As a consequence, NL
is a p edic ed e o and i is a unc ion o he esiduals ob ained
om he linea model:
NL +ph = (e ,e −ph,. . . ,e −n·ph) + ε =ˆ
e +ph, (8)
whe e
ê
is he p edic ed esidual, is he nonlinea unc ion ob ained by he SVR model, nis he size o
he au o eg essi e window, ph is he p edic ion ho izon, and
ε
is he e o e m. Finally, he p edic ion
is ob ained by adding he wo single componen s:
Y +ph =ˆ
L +ph +NL +ph. (9)
O he p oposals in he amewo k o hese kinds o hyb id models ha e eme ged in ecen yea s.
Pa icula ly, inco po a ing he p edic ed alues om a p e ious s ep and he o iginal da a wi h he
esiduals as inpu s o he non-linea o ecas ing model has gained special a en ion. The ime se ies
is conside ed as a unc ion o se e al a iables, such as he o iginal da a, he esiduals, and he
p edic ed alues om he linea model. Some examples o his wo-s ep me hodology a e he wo ks o
Khasei and Bija i [
15
,
21
]. Based on hese las wo ks, he ollowing hyb id app oach 2 was de eloped
(a e applying he SOM me hod) conside ing wo o h ee ypes o a iables as inpu s.
SOM-SARIMA-SVR-2 model (hyb id app oach 2). The ime se ies is he e conside ed a nonlinea
unc ion o he o iginal da a and he esiduals and he p edic ed alues om he i s s ep:
ˆ
Y +ph = (y,e,ˆ
p) + ε +ph, (10)
whe e is he nonlinea unc ion ob ained wi h he SVR model, y
is he o iginal da a, e
is he esidual
ob ained in he SARIMA model, and
ˆ
p
is he p edic ed alue om he SARIMA model a ime .
These a iables a e p esen ed in an au o eg essi e o m, as i is exp essed in he ollowing equa ion:
ˆ
Y +ph = (y ,y −1·ph,y −2·ph . . . ,y −ny·ph,e +ph,e ,e −1·ph,. . . ,e −ne·ph,ˆ
p +ph,ˆ
p ,ˆ
p −1·ph,. . . ,ˆ
p −np·ph) + ε +ph (11)
whe e ne,ny and np ep esen au o eg essi e window sizes o e,yand
ˆ
p
a iables, espec i ely.
Conside ing wo o h ee a iables, and hei au o eg essi e window sizes, a la ge numbe o possible
unc ions can be es ed.
These wo hyb id app oaches a e used in each clus e s o med in he second s ep. Conside ing
wo clus e s, he pa ame e anges used a e also p esen ed in Table 1. No e ha in cases which a iables
eand pa e in ol ed, hei i s alue used as inpu is he p edic ed alue a +ph ime. The same
canno be said o he o iginal da a a iable y, due o he y alue a ime +ph is he inal p edic ion
pu sued. In all ins ances, he SARIMA model is necessa y o ob ain he equi ed da a used as inpu s in
he hi d s ep, oge he wi h he o iginal pa e ns o he con aine se ies.
Appl. Sci. 2020,10, 8326 16 o 17
he combina ion o models and he hyb idiza ion o hem. In his s udy, a combined-hyb id
SOM-SARIMA-SVR o ecas ing model has been p oposed based on a h ee-s ep p ocedu e o p edic
he daily numbe o con aine s passing h ough a Bo de Inspec ion Pos o a ma i ime po .
To educe he complexi y o he p oblem, a clus e ing SOM is i s applied o ob ain smalle
egions wi h simila s a is ical ea u es which may be easie o p edic . A SARIMA model is hen i ed
wi hin each clus e o ob ain p edic ed alues and esiduals o he clus e ed da abase. Finally, a SVR
model is used o o ecas each clus e independen ly using he a iables ob ained om he second s ep
oge he wi h he o iginal da a as inpu s. The combina ion o each clus e esul s in he whole p edic ed
ime se ies. The abo e me hodology in ol es he ad an ages o combining a o ecas ing model wi h a
clus e ing echnique and he s eng hs o he hyb id models in cap u e linea and nonlinea pa e ns.
The p oposed SOM-SARIMA-SVR model has been de eloped and compa ed o o he possible
me hodologies implied in he p ocess (SVR, SOM-SVR and SARIMA-SVR). The esul s showed ha he
SOM-SARIMA-SVR model was he mos compe i i e model, imp o ing he o ecas ing pe o mance
o he es o he models conce ning he p edic ion o he con aine demand, ou pe o ming hese
me hodologies. Pa icula ly, conside ing he SOM-SARIMA-SVR model, wo hyb id app oaches
we e assessed: he classical addi i e app oach, whe e he e o e m (e) o he linea model is he
inpu o he nonlinea model; and he p oposed app oach, whe e he p edic ion is a unc ion o he
o iginal da a (y) and he e o e m (e), and he p edic ion o he linea model (p). Mos accu a e
esul s we e yielded by he second app oach in all cases es ed. In addi ion, he e a e no signi ican
di e ences be ween he p edic ion pe o mance using bo h p edic ion ho izons. The same occu s wi h
he in oduc ion o ce ain a iables as inpu s. These ou comes highligh he obus ness o he model.
This in es iga ion sugges s he eby ha he p oposed hyb id app oach achie es be e ou comes
ge ing highe o ecas ing pe o mance han he classical addi i e hyb id app oach.
To conclude, his s udy is he i s one in using he SOM-SARIMA-SVR model o o ecas he
olume o con aine s passing h ough a BIP o po s in pa icula , and o ime se ies o ecas ing in
gene al wi h p omising esul s. Due o i s abili y o seize he s eng hs o linea and nonlinea models
and he eby o cap u e linea and nonlinea pa e ns, he p oposed model could be applied in o he
ime se ies whe e hese pa e ns a e join ly p esen ed. Pa icula ly, he use o a clus e ing echnique
allows educing he complexi y o he ime se ies, inc easing he accu acy o he inal p edic ion.
Ob iously, some limi a ions a e ound in he model. As wi h any da a-d i en model, when comple ely
di e en inpu s come in o he model, he p edic ion accu acy could wo se. Thus, a e ained model
and a eadjus men o he pa ame e s migh be equi ed o e ime.
Fu u e wo ks will ocus on he applica ion o o he la es non-linea echniques, such as Deep
Lea ning, in o de o assess any imp o emen s in he p edic ion. Knowing he daily con aine demand
in ad ance allows de ec ing wo kload peaks in a po acili y. This gua an ees he co ec planning and
o ganiza ion o a ailable human and ma e ial esou ces. The p oposed me hodology can p o ide an
au oma ic ool o p edic wo kloads a sani a y acili ies a oiding conges ion and delays. The e o e,
i can be used as a decision-making ool by po manage s due o i s capaci y o plan esou ces
in ad ance.
Au ho Con ibu ions: Concep ualiza ion, J.J.R.-A. and I.T.; me hodology, J.J.R.-A.; so wa e, J.J.R.-A., J.A.M.-L.,
D.U.; alida ion, J.J.R.-A., D.U. and J.G.-E.; o mal analysis, J.J.R.-A.; in es iga ion, J.J.R.-A. and J.A.M.-L.; esou ces,
D.U. and J.G.-E.; da a cu a ion, J.J.R.-A. and J.A.M.-L.; w i ing—o iginal d a p epa a ion, J.J.R.-A. and J.A.M.-L.;
w i ing— e iew and edi ing, D.U., J.G.-E. and I.T.; isualiza ion, J.J.R.-A.; supe ision, J.J.R.-A. and I.T.; p ojec
adminis a ion, I.T.; unding acquisi ion, I.T. All au ho s ha e ead and ag eed o he published e sion o
he manusc ip .
Funding:
This esea ch was unded by MICINN (Minis e io de Ciencia e Inno aci
ó
n-Spain), g an numbe
RTI2018-098160-B-I00.
Acknowledgmen s: The da abase was kindly p o ided by he Po o Algeci as Bay Au ho i y.
Con lic s o In e es : The au ho s decla e no con lic o in e es .

Appl. Sci. 2020,10, 8326 17 o 17
Re e ences
1.
Xie, G.; Wang, S.; Zhao, Y.; Lai, K.K. Hyb id app oaches based on LSSVR model o con aine h oughpu
o ecas ing: A compa a i e s udy. Appl. So Compu . 2013,13, 2232–2241. [C ossRe ]
2.
Box, G.E.P.; Jenkins, G.M. Time Se ies Analysis: Fo ecas ing and Con ol; Holden-Day: Oakland, CA, USA, 1976.
3.
Gamba della, L.M.; Bon empi, G.; Tailla d, E.; Romanengo, D.; Raso, G.; Pie ma i, P. Simula ion and
o ecas ing in in e modal con aine e minal. In P oceedings o he 8 h Eu opean Simula ion Symposium;
SCS In e na ional: Genoa, I aly, 1976; pp. 626–630.
4.
Babcock, M.W.; Lu, X. Fo ecas ing inland wa e way g ain a ic. T ansp. Res. Pa E Logis . T ansp. Re .
2002
,
38, 65–74. [C ossRe ]
5.
Klein, A. Fo ecas ing he An we p ma i ime a ic lows using ans o ma ions and in e en ion models.
J. Fo ecas . 1998,15, 395–412. [C ossRe ]
6.
Peng, W.-Y.; Chu, C.-W. A compa ison o uni a ia e me hods o o ecas ing con aine h oughpu olumes.
Ma h. Compu . Model. 2009,50, 1045–1057. [C ossRe ]
7.
Zhang, G.P. Time se ies o ecas ing using a hyb id ARIMA and neu al ne wo k model. Neu ocompu ing
2003
,
50, 159–175. [C ossRe ]
8. Kohonen, T. Sel O ganising Maps; Sp inge : Be lin, Ge many, 1995.
9.
De Bod , E.; Co ell, M.; Le emy, P.; Ve leysen, M. On he use o sel -o ganizing maps o accele a e ec o
quan iza ion. Neu ocompu ing 2004,56, 187–203. [C ossRe ]
10.
Vlahogianni, E.I.; Ka la is, M.G.; Golias, J.C. Tempo al E olu ion o Sho -Te m U ban T a ic Flow:
A Nonlinea Dynamics App oach. Compu . Ci . In as uc . Eng. 2008,23, 536–548. [C ossRe ]
11.
Chen, H.; G an -Mulle , S.; Mussone, L.; Mon gome y, F. A s udy o hyb id neu al ne wo k app oaches and
he e ec s o missing da a on a ic o ecas ing. Neu al Compu . Appl. 2001,10, 277–286. [C ossRe ]
12.
Ismail, S.; Shab i, A.; Samsudin, R. A hyb id model o sel -o ganizing maps (SOM) and leas squa e suppo
ec o machine (LSSVM) o ime-se ies o ecas ing. Expe Sys . Appl. 2011,38, 10574–10578. [C ossRe ]
13.
Ma, Z.; Xing, J.; Mesbah, M.; Fe ei a, L. P edic ing sho - e m bus passenge demand using a pa e n hyb id
app oach. T ansp. Res. Pa C Eme g. Technol. 2014,39, 148–163. [C ossRe ]
14.
Tseng, F.-M.; Yu, H.-C.; Tzeng, G.-H. Combining neu al ne wo k model wi h seasonal ime se ies ARIMA
model. Technol. Fo ecas . Soc. Chang. 2002,69, 71–87. [C ossRe ]
15.
Khashei, M.; Bija i, M. A no el hyb idiza ion o a i icial neu al ne wo ks and ARIMA models o ime se ies
o ecas ing. Appl. So Compu . 2011,11, 2664–2675. [C ossRe ]
16.
Wang, J.-J.; Wang, J.-Z.; Zhang, Z.-G.; Guo, S.-P. S ock index o ecas ing based on a hyb id model. Omega
2012,40, 758–766. [C ossRe ]
17.
Chen, K.-Y.; Wang, C.-H. A hyb id SARIMA and suppo ec o machines in o ecas ing he p oduc ion
alues o he machine y indus y in Taiwan. Expe Sys . Appl. 2007,32, 254–264. [C ossRe ]
18.
Zhu, B.; Wei, Y. Ca bon p ice o ecas ing wi h a no el hyb id ARIMA and leas squa es suppo ec o
machines me hodology. Omega 2013,41, 517–524. [C ossRe ]
19.
Kohonen, T. Sel -o ganized o ma ion o opologically co ec ea u e maps. Biol. Cybe n.
1982
,43, 59–69.
[C ossRe ]
20. Vapnik, V.N. S a is ical Lea ning Theo y; John Wiley and Sons: New Yo k, NY, USA, 1998.
21.
Khashei, M.; Bija i, M. An a i icial neu al ne wo k (p, d, q) model o imese ies o ecas ing. Expe Sys . Appl.
2010,37, 479–489. [C ossRe ]
Publishe ’s No e:
MDPI s ays neu al wi h ega d o ju isdic ional claims in published maps and ins i u ional
a ilia ions.
©
2020 by he au ho s. Licensee MDPI, Basel, Swi ze land. This a icle is an open access
a icle dis ibu ed unde he e ms and condi ions o he C ea i e Commons A ibu ion
(CC BY) license (h p://c ea i ecommons.o g/licenses/by/4.0/).