Stochastic Fluid Flow Model Analysis
Markus Fiedler 1
andHolgerVoos 2
1
UniversityofKarlskrona/Ronneby
Dep.ofTelecommunicationsandSignalProcessing (ITS)
S-37179Karlskrona,Sweden
phone:+46-708-537339, fax:+46-455-385667 markus.fiedler@its.hk-r.se
2
UniversityofKaiserslautern,InstituteofProcessAutomation
PO-Box3049,D-67653 Kaiserslautern,Germany
phone:+49-631-2054457,fax:+49-631-2054462,voos@eit.uni-kl.de
Abstract. Thestochastic uid owmodel(SFF)isoneoftheleading
modelsinperformanceevaluationfortele-anddatacommunicationsys-
tems,especially infastpacket-switching networksand ATM.However,
thenumerical analysisof the SFFis widelyconsideredto be unstable.
Inthis paper,someinvestigationsandresults arepresented concerning
thenumerical stability of the SFFanalysis also forlarge systemswith
nitebu er.Weidentifythemainsourceofthenumericalproblemsand
give hintshowtocircumventthem.Theusefulnessof di erent solution
methodsare comparedand the mostrobust methodsfor systemswith
largenumbersofsourcesandlargebu ersizesareidenti ed.
1 Introduction
Onepossiblemodelforperformanceevaluationincommunicationssystemsisthe
socalledstochastic uid owmodel(SFF). Itiswidelyused whendealingwith
fastpacket-switchedandATMnetworksbutalsoallowsperformanceanalysisin
anykindoftele-anddatacommunicationssystemincludingmobilecommunica-
tions.Firstpioneering work consideringtheSFFcan befoundin [8].In[1],an
eleganttreatmentofa nite numberofhomogeneouson-o uidsourceswhose
streamswereconcentratedbyamultiplexerwithin nitebu erispresented.In
1984,Kostenpresentedanexpansionofthat modeltoheterogeneoustraÆc[9],
andin1988,Tuckerpublished theformalwayofhowtotreat nitemultiplexer
bu ers [17].The main advantage of the SFF is the computational complexity
whichisindependentofthebu ersize[12].
However, the analysis of the SFF is widely considered to be numerically
unstable. Indeed, the solutioncontainscomponents that are exponentiallyin-
creasingwiththebu ercontentinthecaseoflimitedbu ersize[8],[1],[17].For
that reason andthe fact that no closed-formsolutionsexist,most work in the
1990'sdealswithbu ersofin nitesizeorapproximationsofthesolutionbased
upon that assumption. In 1991, Nagarajan et al [13] reported results for 130
unable to obtain resultsfor large-sizeproblems due to numericalinstability of
theirsolutionmethod.InspiteofthesediÆculties,veryfewmaterialexist that
discussesthesenumericalproblemsindetail.
Therefore,thispaperclosesthatgapbetweentheoryandnumericalpractice
fortheclassicalspectralmethodthatisbasedoncomputinganeigensystem(spec-
traldecomposition andondeterminingthecoeÆcientsassociatedtothespectral
components.Wearegoingtodemonstratehowstablethespectralmethodactu-
allymightperformevenforlargesystemswith nitebu ersifsomequitesimple
rules are taken care of. In section 2, we describe the model and the notation
that is used throughoutthis paper.Section 3dealswith themain stepsof the
stochastic uid ow analysis. Section 4 introduces the numerical methods to
solve the system of linear equations, which is necessaryto adapt the solution
ofthesystemofdi erentialequationstotheboundaryconditions.Furthermore,
somespecialimplementation issuesare proposed.Section 5presentsnumerical
resultsand discusses thefeasibility of di erentnumerical methods forsystems
with homogeneousand heterogeneous sources. Section 6summarizes themain
resultsandgivesanoutlookonfuturework.
2 The stochastic uid ow model
TheSFFunderinvestigationistheclassicalmodelwhichcanbefoundinmany
publications.WeassumeN on-o uidsourcesthatalternaterandomlybetween
an on state with peak cell rate h and an o state, both of exponentiallydis-
tributed duration. Thus, a discrete-valued rate process S is formed with n
S
statess
i
and ow intensities R 2fr
i
g;i 2f0;:::;n
S
1g. Inaddition, we as-
sume a uid bu er of nite size K and an outlet (= server) with capacity C.
ThedynamicoftherateprocessSisdescribedbyanirreducibleMarkovchain,
representedbythein nitesimalgeneratormatrixM,whichalsodeterminesthe
state probabilities
i
. Each statehas a drift valued
i
= r
i
C, depending on
whichthestatesareclassi edin
{ over-loadstatesS o
=fi
d
i
>0g;
{ equilibriumstatesS e
=fi
d
i
=0g;
{ under-loadstatesS u
=fi
d
i
<0g:
ThesevaluesarecollectedinthesocalleddriftmatrixD=diag [d
i
].Thematrix
R = diag [r
i
] is called rate matrix. We denote the aggregatestate space with
S =S o
[S e
[S u
.Thenon-o sourcesmightbehomogeneous,i.e.allhavethe
sameparameters,orheterogeneous.Inthelattercase,weformn
G
groups,each
containinghomogeneoussources.
3 The uid ow analysis
Thefollowingrepresentsasummary ofthe uid owanalysiswhichdealswith
the mostessentialpoints in the numerical context. Further information might
stationarydistributionfunctionF(x)withelementsF
i
(x)=PrfXx^state=
igisgovernedbythesystemof di erentialequations
D d
dx
F(x)=MF(x): (1)
From(1), aneigenvalue-eigenvectorproblem
z
q D'
q
=M'
q
(2)
with q 2 f0;:::;n
S
1g is obtained. Theeigenvalues z
q
and eigenvectors'
q
appearinthesocalled spectralcomponentsinthesolutionof(1),
F(x)= X
q2S a
q (K)'
q exp(z
q
x): (3)
Forhomogeneouson-o sources,eigenvaluesandeigenvectorsaregiveninclosed
form [1]. We normalizethe eigenvectorsin awaythat the sumof all elements
equalsone.If thesourcesare heterogeneous,theset ofeigenvaluesfz
q
ghas to
bedetermined numericallyfrom theinverseeigenvalueproblem
q (z
q )'
q
=
R 1
z
q M
'
q
with
q (z
q )=
n
G
X
j=1
(j)
q (z
q
)=C: (4)
Forgroupsofon-o sources,thefunctions (j)
q (z
q
)arealsogivenin closedform
anddependonthenumberofsourcesthatareoninstateqandbelongtogroupj.
Once z
q
isdetermined,the correspondingeigenvectorisacompositionof parts
that are determined for homogeneous groups. This composition is done using
Kroneckeralgebra.Moredetails mightbefoundin[15],[10].
The coeÆcients a
q
(K) in (3) are necessary to adjust the solution to the
boundaryconditions,whichdepend onthebu ersizeK:
F
i
(K)=
i
; i2S u
; F
i
(0)=0; i2S o
: (5)
Insertion of (5) into (3) leads to a system of linear equations that has to be
solvedinorder toobtainthecoeÆcients:
X
q2S u
[S o
a
q (K)'
qi exp(z
q K)=
i
; i2S u
X
q2S u
[S o
a
q (K)'
qi
=0; i2S o
(6)
Originally,casesd
i
=0hadtobeexcludedbychoosinganappropriatevalue
ofC[1].However,weareabletodealwiththoseequilibriumstatesinthesame
spondingcoeÆcientalsoescapesinthelimitwhilethecorrespondingeigenvector
approachesaunityvectorin q-direction:
lim
d
q
!0 a
q
(K)=0; lim
d
q
!0 '
q
=e
q
(7)
This means that in the limitd
q
! 0,there is no couplingbetween equation q
andtheotherequations. Thus, wearefreetoignorestateswithvanishingdrift
duringthesolutionprocedure.
3.1 Finitebu er case
The mostly ill-conditioned system of linear equations (6) of size (dimS u
+
dimS o
) n
S
has to be solved numerically. Di erent methods to do this will
bepresented andevaluated in sections4and 5.Assoonasthe coeÆcientsare
determined,theprobabilitythat thebu eris fullinstateiisfoundtobe
u
i
(K)=
i lim
b!K F
i (b)=
i X
q2S a
q (K)'
q exp(z
q
K); (8)
see[17].Withthis, thelossprobabilitycanbeexpressedas
P
L (K)=
P
i2S o
u
i (K)d
i
E[R ]
: (9)
3.2 Bu er-less uid ow model
Thesocalled bu er-less uid owmodel isobtainedbypassingthebu er size
to the limit K ! 0. This case may be assumed if the bu er is much smaller
thanthemeanburstlengthandthus(almost)loosesitsin uenceonburstlevel.
Inthis case,the uid owanalysis describedbeforedoesn'tneedto becarried
out anymore which allows the computation of systems with very large state
spaces. Losshappens assoon aspositivedriftoccurs, so that theprobabilities
u
i
inthelossprobabilityformula(9)canbereplacedbythecorrespondingstate
probabilities
i :
P
L (0)=
P
i2S o
i d
i
E[R ]
: (10)
4 Numerical methods
In thissection we present aselectionof methods that havebeenused to solve
the system of linear equations (6) which is mostly ill-conditioned, not sparse,
not symmetrical and not positive de nite. In addition, the systems are large
(dependingonthesizeofthestatespaceS)andassociatedwithnumericalinput
errorsthat stemfromthenumericaldeterminationoftheeigensystem.Ouraim
is to checkhowthefollowingmethods work in spite ofthese diÆculties. Since
theclassi cationofthemethods isnotuniqueintheliterature,werefertothe
With direct methods, we describe \numerical methods that computesolutions
of mathematical problems in a xed number of operations" [16]. We focus on
methods that transform asystem of linearequations Ax = binto asystem
A 0
x=b 0
during aso-calledreductionphase. Herein,A 0
representsanupper
triangularmatrix.Then, thecomponentsof thesolutionvectorxarefoundby
backsubstitution.See[5],[16],[4]fordetailsonthefollowingmethods:
1. Gaussian elimination.The reductionis achieved by aspecial succession of
divisionsandsubtractionsofrows.Assume thatin stepithe rsti 1rows
havealreadybeentreated.Theelementsofrowiareobtainedby
a (i)
k l
=a (i 1)
k l a
(i 1)
k i
a (i 1)
ii a
(i 1)
il
k>i;l=1:::n: (11)
Theelementsa (i 1)
ii
arecalledpivots.
2. Gaussianelimination withpartialpivoting. Fromanumericalpointofview,
theelementsonthediagonalmightnotrepresentoptimal pivots.It canbe
shown that theabsolute valueof thepivot should be as large aspossible.
Hence,if apivot withlarger absolute value canbeobtained from another
rowin thesamecolumn,thenthecorrespondingrowsareinterchanged.
3. Gaussianeliminationwithfullpivoting.Here,anoptimalpivotissearchedin
rowsandcolumns.However,interchangingofcolumnsleadstointerchanged
elementsinthesolutionvector(asystemA 0
x 0
=b 0
appears), whichhave
tobere-changedafter backsubstitution.
4. Givens rotations. Thereductionis performed bymultiplying thesystemof
equationswithelementaryrotationmatrices,herebycancellingouttheele-
mentsinthelowertriangularpart.In[16],thismethodcanbefoundamong
projection methods.
5. Householdertransformation.Thismethodalsousestransformationmatrices
forthereductionphase.Fordetails,see[11].
Themethods3to5arestablewithregardtothealgorithmerror.
4.2 Iteration methods
Iteration methods try to approximate the solution by carrying out iterations,
consistingofacoupleofoperations,aslongastheapproximationx
hasnotyet
convergedtoadesiredextent.Startingfromx
(0)
=0,theresultofaniteration
servesasbasisfor thenextiteration, i.e. x
(i+1)
=f(x
(i)
).Toget apositive
de nite matrix instead of A that guarantees convergence [5], we rst apply
a so called Gaussian transformation by multiplying our system of equations
with the transposed matrix A T
from the left side, i.e. we treat the system
A T
Ax=A T
b.Seeagain[5],[16],[4]fordetailsonthefollowingmethods,
whichweuseinamoredirectway,aswe xthe(maximal)numberofiterations
x
(k +1)
i
=x
(k )
i +
"
a
ii
0
@
b
i i 1
X
j=1 a
ij x
(k +1)
j
n
X
j=i a
ij x
(k )
j 1
A
(12)
with"=1.
2. Incomplete relaxation. The speed of the convergence of the Gauss-Seidel
method, which tries to minimize (= relax) the defect vector b Ax
completely, mightbeimprovedbyover-relaxationwith ">1.However,the
relaxation factor " has to be chosen carefully; it depends heavily on the
systemofequationsunderstudy.
3. Jacobi rotation. This method uses the same kind of rotation matrices as
Givens rotationsto \improve" thematrix A T
A by makingit diagonally
dominant.Thegradeofthisdominanceisspeci edbytherotationlimitand
therotationstopsassoonasthedesiredgradeisreached.Inthiswork,the
algorithm stops in any caseafter 1000iterations. Finally, adirect method
deliversthesolution.
All thesemethodsaresensitivetoroundingerrors.
4.3 Implementationissues
A proper numerical evaluation of the eigensystem has to be ensured in order
to keep the input error for the solution of the system of linear equations (6)
as small aspossible.Therefore, the eigenvaluesshould be obtained with great
precisionwhichisnoproblemforon-o sources,evenifanumericalsearchbased
on(4) wasperformed.Butalso theeigenvectors'
q
should bedetermined with
care even if the orders of magnitude of certain elements are extremely small.
Section5.3demonstrateswhathappensiftheprecisionoftheeigenvectorcalcu-
lationdrops.Notethatduetolim
di!0 z
i
=1,stateswithsmalldriftsmight
causenumericalunder-/over owin (3).But(7) showsthat thecontribution of
such components becomes arbitrarily small. Hence,they can be extracted be-
foresolvingthesystemoflinearequations(6). As oating-pointvariables with
doubleprecisionrange from about10 323
to about10 308
, wedeleted such
equations whose presence would probably leadto over owin column q due to
exp(z
q
K)>10 300
ortozerosduetoexp(z
q
K)<10 300
.Afterthisreduction,a
numberofn 0
S
equationsisleft.
As public and commercial tools for numerical mathematics mostly do not
allow theuser to look inside their implementation, wewroteowncodefor the
solution methods in C++ with great care to avoid bad numerical surprises.
Toavoidexcessivememoryconsumptionand executingtimes(c.f. section5.3),
oating-point variables with double precision (8 Bytes) were used instead of
Inthissection, wepresentthequalityofnumericalresultsforlossprobabilities
inquitelargesystemswithhomogeneousandheterogeneoussources,whichhave
beenobtainedoncomputersofthetypesSunSparcStation10and20(architec-
turesun4m).Resultsthatwereobviouslywrong,e.g.negativevalues,havebeen
markedwith\|".
Inrelationtothemeanburstsizesoftheon-o sourcesunderconsideration,
bu ersizesareclassi edasfollows:
1. XS:Verysmall bu erthat hasnoin uence onburstlevel.
2. S:Smallbu erof 0.1timestheminimalmeanburstsize.
3. M:Medium-sizedbu eroftheminimalmeanburstsize.
4. L:Largebu erof10timestheminimalmeanburstsize.
5. XL:Verylargebu erof100timestheminimalmeanburstsize.
Incase1,numericalresultsmightbecomparedwithverystableresultsprovided
bythebu er-less uid owmodel.Therefore,resultsdeliveredbythe uid ow
analysisareregardedtobegoodenoughifthecorrespondingrelativeerrordoes
not surpass0.1 %.However, in the cases2 to 5, uid ow simulations had to
beused inorderto obtainreferencevalues.Here,theerrortolerance wasset to
the95% con denceinterval,whosesizewasmostlylessthan4%of thecorre-
spondingaveragevalue.Amethodforsolvingthesystemoflinearequations(6)
isconsideredasapplicableforoneofthefollowingsystemsifthecorresponding
errortoleranceiskept.
5.1 Homogeneous system
The homogeneous system under study consists of N 2 f50;100;:::;600g quite
bursty on-o sources with amean-to-peak bit rate ratioof 0.1. These sources
are muchmore bursty and thus morecritical than thevoicemodelin [7]. The
capacityofthemultiplexerwaschoseninawaythattheloadequals0.8;thesize
ofthestatespaceisn
S
=N+1.
XS bu er:Table1comparesresultsobtainedbyusingonedirectandtwoiter-
ativemethodswithreferenceresultsobtainedfromthebu er-lessanalysis.The
Table1.Lossprobabilitiesobtainedwithsomesolutionmethods,givenhomogeneous
sourcesandaverysmallbu er.
N Fullpivot. Gauss-Seidel Incompl.relax.(factor)Jacobi/Bu er-less
100 3:738910 2
3:623510 2
3:695510 2
(1.9) 3:738910 2
200 1:376810 2
1:353210 2
1:371110 2
(1.9) 1:376810 2
300 6:687210 3
6:606910 3
6:687510 3
(1.9) 6:687210 3
400 3:485210 3
3:459010 3
3:485110 3
(1.4) 3:485210 3
500 1:987510 3
1:979910 3
1:989410 3
(1.2) 1:987510 3
600 | 1:147510 3
1:147510 3
(1.0) 1:148010 3
Householder method already fails if the number of sources exceeds 300. The
Gauss-Seidel method produces better results the more the size of the system
increases.Theincompleterelaxationmethodleadsto furtherimprovement,but
only ifthe relaxationfactor hasbeenoptimized before. As there are reference
valuesavailable,wechosethebestresultforeachN basedonrelaxationfactors
2f1:0;1:1;:::;1:9g.TheJacobimethod witharotationlimitof0.9reproduced
the reference values exactly for any number of sources, no matter which di-
rectmethod is usedat itsend. For anumberof600sources,onlytheiteration
methodsarestillapplicable.Table2summarizesthecasesinwhichthedi erent
methodsproducedresultswithinthede nederrortoleranceof0.1%;thesecases
aremarkedby\?".
Table 2.Usefulness of the solutionmethods given homogeneoussources and avery
smallbu er.
Gausswithoutpivoting ? ? ? ? ? ? ? ? ?
Gausswithpartialpivoting ? ? ? ? ? ? ? ? ? ? ?
Gausswithfull pivoting ? ? ? ? ? ? ? ? ? ? ?
Givensrotations ? ? ? ? ? ? ? ? ? ? ?
Householder ? ? ? ? ? ?
Gauss-Seidel ?
Incompleterelaxation ? ? ? ? ? ? ?
Jacobirotation ? ? ? ? ? ? ? ? ? ? ? ?
Numberofsources 100 200 300 400 500 600
M bu er: Table3showssomeresultsobtainedwithdirectmethods andcom-
pares them withsimulationresults. For 550sources, the Gaussian elimination
Table 3.Lossprobabilitiesobtainedwithsomesolutionmethodsgivenhomogeneous
sourcesandamedium-sizedbu er.
N Gauss Fullpivot. Givensrot. Simulation
300 3:554310 3
3:554410 3
3:554410 3
(3:57010:1202)10 3
350 2:504510 3
2:504410 3
2:504410 3
(2:49810:0864)10 3
400 1:799210 3
1:798010 3
1:798010 3
(1:78190:0553)10 3
450 1:313010 3
1:313310 3
1:313410 3
(1:30340:0550)10 3
500 9:683810 4
9:702610 4
9:730210 4
(9:71760:3259)10 4
550 7:616410 4
7:242810 4
7:223010 4
(7:06910:2103)10 4
600 | | 5:456010
4
(5:27170:2138)10 4
withoutpivoting deviatesobviously,while Givensrotationsistheonlymethod
that manages 600 sources with acceptable precision. As in the previous case,
theHouseholdermethodisnotabletotreatmorethan300sources.Amongthe
iterativemethods, most of the Gauss-Seidelresults(not shown) lie within the
corresponding con dence interval. Due to the lack of reference values, i.e. the
possibility to optimizetherelaxationfactor,the incompleterelaxationmethod
are summarized in Table4, where \?" stands for aresult that lies within the
correspondingcon denceinterval.
Table4.Usefulnessofthesolutionmethodsgivenhomogeneoussourcesandamedium-
sizedbu er.
Gausswithoutpivoting ? ? ? ? ? ? ? ? ? ?
Gausswithpartialpivoting ? ? ? ? ? ? ? ? ? ? ?
Gausswithfull pivoting ? ? ? ? ? ? ? ? ? ? ?
Givensrotations ? ? ? ? ? ? ? ? ? ? ? ?
Householder ? ? ? ? ? ?
Gauss-Seidel ? ? ? ? ? ? ? ?
Jacobirotation ? ? ? ? ? ? ? ?
Numberofsources 100 200 300 400 500 600
L bu er: ForL bu er, thebehaviorof the direct methods is quite similar to
thatintheMbu ercase.Wecon neourselvestothesummarythatispresented
inTable5.TheGauss-SeidelmethodshowslargerdeviationsasintheMbu er
Table 5. Usefulness ofthesolutionmethodsgivenhomogeneous sourcesand alarge
bu er.
Gausswithoutpivoting ? ? ? ? ? ? ? ? ?
Gausswithpartialpivoting ? ? ? ? ? ? ? ? ?
Gausswithfull pivoting ? ? ? ? ? ? ? ? ? ? ?
Givensrotations ? ? ? ? ? ? ? ? ? ? ? ?
Householder ? ? ? ? ? ?
Gauss-Seidel ?
Jacobirotation
Numberofsources 100 200 300 400 500 600
case;merelytheresultfor 300sourceslieswithin thecorrespondingcon dence
interval,while themethodfailsfor morethan350sources.TheJacobirotation
wasnotabletoproduceanymeaningfulresult.ForLbu er,thedirectmethods
seemtobemuchbetterapplicablethantheconsiderediterationmethods.
5.2 Heterogeneoussystem
Here,weassumetwogroupsof30on-o sources,each withthe samemean bit
rateanddi erentparametersasfollows:
1. Mean-to-peak bit rate ratio of 0.5 (less bursty), ratio of bu er and mean
burstsizeK =hb
1
2f0;1;10;100g;
2. Mean-to-peakbit rate ratioof 0.1 (morebursty), ratioof bu erand mean
burstsizeK =hb
2
2f0;1;10;100g.
For XS bu er, each direct method is able to compute the corresponding
lessmodel. ThisbehaviormightbeinterpretedasasignalfortheGauss-Seidel
method's sensitivitytoinputerrors.
Table 7 shows results for the Gaussian methods for di erent bu er sizes.
The result for XLbu er is too small to be compared with simulation results.
A comparison with the upper bound for the loss probability that is based on
thedominanteigenvalue[6],[2]excludestheresultsobtainedwiththeGaussian
methodwithoutandwithpartialpivoting;thesameisvalidforGivensrotations
resultsthat arenotshownexplicitly. As thisupper bound isknownto overes-
timate thereal resultby severalorders of magnitude,the resultobtainedwith
the Gaussian method with full pivoting seemsreasonable. This method seems
tobebestamongthedirectmethodswhendealingwithsmallprobabilities.On
theotherhand,HouseholderandGauss-Seidelfail incasesoflargebu ersizes.
Table8showsasummary.
5.3 The role ofthe eigenvectors
In this section, we show to which extent the accuracy of the eigenvectorsdi-
rectlyin uencesthenumericalstabilityof thesolutionmethods.Here,onlythe
homogeneous systems need to be investigated (in the heterogeneous case, the
eigenvectorsareobtainedasacomposition ofparts thataredeterminedforho-
mogeneousgroups).Sinceinthis homogeneouscaseclosedformuladescriptions
of the eigenvaluesexist, theycan be computed with great precision. However,
duringthecalculationoftheeigenvectors,hugedi erencesintheordersofmagni-
tudeoccur:thereareproductsthatcontaintermscausingnumericalover owas
wellasnumericalunder owwhiletheresultingproductwouldlieinthetractable
rangeconcerningtheorderofmagnitude.Inoursystem,theuseof oating-point
numberswithdoubleprecision(8Bytes)within thatcriticalproductrestricted
thenumberofsourcesto100.
A rst improvementwas reached by a separate treatment of base and ex-
ponentofeigenvectorelementsintwo oating-pointnumberswithdoublepre-
cision, i.e. a kindof logarithmiccalculation within the critical product. Using
this trick, thenumberof sourcescould beraisedto 250.Theprice tobe paid:
Onaverage,thecalculationoftheeigenvectorstookmorethan10timesaslong.
However,theresultsforupto600sourcesthatwerereportedintheprevioussec-
tion havebeen obtainedusing long double oating-point numbers (16Byte).
Compared to the useof the standarddoubledata type,the speed went down
Table6.LossprobabilitiesobtainedwiththeGauss-Seidelmethodforheterogeneous
sourceswithdi erentmeanburstsizesandverysmallbu er.
hb
1 hb
2
Gauss-Seidel Reference
10 10 2:428810 3
2:437610 3
10100 2:426010 3
2:437610 3
100 10 2:448910 3
2:437610 3
3 3
neoussourcesandMtoXLbu ers.
K
hb
1 K
hb
2
Nopivoting Partialpivot. Fullpivot. Sim./Approx.
1 1 4:264810 4
4:264810 4
4:264810 4
(4:21320:1059)10 4
10 10 1:141510 7
1:141510 7
1:141510 7
(1:12400:0819)10 7
100100 1:077810 18
7:500610 19
1:512010 37
Approx.:1:299210 33
Table 8. Usefulnessofthesolutionmethodsfor di erent bu ersizesgivenheteroge-
neoussources.
Gausswithoutpivoting ? ? ?
Gausswithpartialpivoting ? ? ?
Gausswithfullpivoting ? ? ? ?
Givensrotations ? ? ?
Householder ? ?
Gauss-Seidel ? ?
Bu ersize: XS M L XL
byafactorofmorethan100.Ingeneral,thee ectofinputerrorsonnumerical
stabilityseemto growasthebu erbecomeslarger.Due to(3), thisisnotsur-
prising:Theeigenvectorsaremultipliedbyexponentialfunctionscontainingthe
bu ersize.Especiallyiterationmethodsseemto su erfromthisproblem; even
the powerful Jacobirotation method fails. The Gaussian elimination with full
pivoting seemstobethebest possiblemethodto getalongwiththisproblem.
6 Conclusions
Inthispaper,wepresentedsomerecipesonhowtostabilizethestochastic uid
owanalysis numericallyif nite bu ers and theclassicalspectralmethod are
used. The most essential step is to determine the eigensystem as precisely as
possible.Thisshouldbedonetokeeptheinputerrorforthenextsteplow,which
isthesolutionofasystemoflinearequationsinordertoadaptthesolutionofthe
systemofdi erentialequationstotheboundaryconditions.Beforesolvingthat
system,thoseequationsthatbelongtostateswithverysmalldriftvaluesshould
beextracted.Thesestateshavenosigni cantin uenceonthesolutionbut may
cause numerical under- or over ow. To solve the system of linear equations,
the solution method has to be chosen very carefully. Considering the list of
possiblesolutionmethodsthatweretestedinthiswork,theuseoftheGaussian
algorithmwithcompletepivotingledtothebestresults,followedbytheGivens
rotationsmethod.Althoughthesummarytablespresentedinsection5arevalid
only for the corresponding parameter settings, they re ect the experience we
also made with other systems. Altogether, the proposed measures lead to a
numericalperformanceof thestochastic uid owmodelanalysis thatis much
betterthanitsreputation.Futureworkonmumericalproblemsinthe uid ow
modelcontextshouldalsoconsiderthematrix-analyticsolutionbyRamaswami
[1] D. Anick, D. Mitra and M. M. Sondhi. Stochastic theory of a data-
handlingsystemwithmultiplesources.The BellSystem TechnicalJournal,
61(8):1871{1894(1982).
[2] A.I.ElwalidandD.Mitra.E ectivebandwitdhofgeneralMarkoviantraÆc
sourcesandadmissioncontrolofhighspeednetworks.IEEE/ACMTransac-
tionsonNetworking,1(3):329{343(1993).
[3] A.I. Elwalid and D. Mitra. Statistical multiplexing with loss priorities in
rate-basedcongestioncontrolofhigh-speednetworks.IEEETransactionson
Communications,42(12):2989{3002(1994).
[4] Encyclopaediaofmathematics.Kluwer,1995,ISBN1-55608-010-7.
[5] D.K.FaddeevandV.N.Faddeeva.Computationalmethodsoflinearalgebra.
Freeman,1963.
[6] R.Guerin,H.AhmadiandM.Naghshineh.Equivalentcapacityanditsap-
plicationtobandwidthallocationinhigh-speednetworks.IEEEJournalon
SelectedAreasinCommunications,9(7):968{981 (1991).
[7] H. He es and D.M. Lucantoni. A Markov modulated characterization of
packetized voiceand datatraÆc and related statistical multiplexerperfor-
mance. IEEE Journal on Selected Areas in Communications, 4(6):856{867
(1986).
[8] L.Kosten.Stochastictheoryofamulti{entrybu er.DelftProgress Report,
SeriesF,1:10{18(1974).
[9] L.Kosten.Stochastictheoryofdatahandlingsystemswithgroupsofmultiple
sources. Proc. of IFIP WG 7.3 { TC 6 International Symposium on the
Performance of Computer Communication Systems, Zurich,1984, pp 321{
331.
[10] K.P. Kontovasilis and N.M. Mitrou. Bursty traÆc modelling and eÆcient
analysisalgorithmsvia uid- owmodelsfor ATMIBCN.AnnalsofOpera-
tionsResearch,49:279{323 (1994).
[11] P. Lancasterand M. Tismentzky.The Theory of Matrices, Second edition
withApplications.New York:AcademicPress, 1985.
[12] D.Mitra.Stochastictheoryofa uid owmodelofproducersandconsumers
coupledbyabu er.AdvancesinAppliedProbability,20:646{676 (1988).
[13] R.Nagarajan, J.F.Kurose andD. Towsley.Approximationstechniquesfor
computing packet loss in nite-bu ered voice multiplexers. IEEE JSAC,
9(3):368{377(1991).
[14] V. Ramaswami. Matrix analytic methods for stochastic uid ows. Proc.
ITC-16,Edinburgh,1999,pp1019{1030.
[15] T. Stern and A. Elwalid. Analysis of separable markov-modulated rate
models for information-handling systems.Advances in Applied Probability,
23:105{139(1991).
[16] W.J. Stewart. Introduction to the Numerical Solution of Markov Chains.
PrincetonUniversityPress, 1994,ISBN0-691-03699-3.
[17] R.Tucker.Accuratemethodforanalysisofapacket-speechmultiplexerwith
limiteddelay.IEEETransactionsonCommunications,36(4):479{483(1988).
[18] T. Yang and D.H.K. Tsang. A novel approach to estimating the cell loss
probabilityinanATMmultiplexerloadedwithhomogeneouson-o sources.
IEEETransactionsonCommunications,43(1):117{126(1995).
