LINE Solver
MATLAB API documentation
Loading...
Searching...
No Matches
FluidPrioQueue.m
1% Ret = FluidPrioQueue(Q, Rin, d, ...)
2%
3% Returns various performane measures of a continuous time
4% fluid priority queue, see [1]_.
5%
6% Parameters
7% ----------
8% Q : matrix of shape (N,N)
9% The generator of the Markov chain modulating the arrival process
10% Rin : matrix of shape (K,N)
11% The matrix defining the input fluid rates in various states of the
12% background process of all fluid types
13% d : real number
14% The state independent fluid service rate
15% further parameters :
16% The rest of the function parameters specify the options
17% and the performance measures to be computed.
18%
19% The supported performance measures and options in this
20% function are:
21%
22% +----------------+--------------------+----------------------------------------+
23% | Parameter name | Input parameters | Output |
24% +================+====================+========================================+
25% | "flMoms" | Number of moments | The moments of the fluid level |
26% +----------------+--------------------+----------------------------------------+
27% | "flDistr" | A vector of points | The fluid level distribution at the |
28% | | | requested points |
29% +----------------+--------------------+----------------------------------------+
30% | "stMoms" | Number of moments | The sojourn time moments of fluid |
31% | | | drops |
32% +----------------+--------------------+----------------------------------------+
33% | "stDistr" | A vector of points | The sojourn time distribution at the |
34% | | | requested points (cummulative, cdf) |
35% +----------------+--------------------+----------------------------------------+
36% | "prec" | The precision | Numerical precision used as a stopping |
37% | | | condition when solving the Riccati and |
38% | | | the matrix-quadratic equations |
39% +----------------+--------------------+----------------------------------------+
40% | "erlMaxOrder" | Integer number | The maximal Erlang order used in the |
41% | | | erlangization procedure. The default |
42% | | | value is 200. |
43% +----------------+--------------------+----------------------------------------+
44% | "classes" | Vector of integers | Only the performance measures |
45% | | | belonging to these classes are |
46% | | | returned. If not given, all classes |
47% | | | are analyzed. |
48% +----------------+--------------------+----------------------------------------+
49%
50% Returns
51% -------
52% Ret : list of the performance measures
53% Each entry of the list corresponds to a performance
54% measure requested. Each entry is a matrix, where the
55% columns belong to the various job types.
56% If there is just a single item,
57% then it is not put into a list.
58%
59% References
60% ----------
61% .. [1] G. Horvath, "Efficient analysis of the MMAP[K]/PH[K]/1
62% priority queue", European Journal of Operational
63% Research, 246(1), 128-139, 2015.
64
65function varargout = FluidPrioQueue(Q, R, d, varargin)
66
67 K = size(R,1);
68
69 % parse options
70 erlMaxOrder = 200;
71 precision = 1e-14;
72 classes = 1:K;
73 eaten = [];
74 for i=1:length(varargin)
75 if strcmp(varargin{i},'erlMaxOrder')
76 erlMaxOrder = varargin{i+1};
77 eaten = [eaten, i, i+1];
78 elseif strcmp(varargin{i},'prec')
79 precision = varargin{i+1};
80 eaten = [eaten, i, i+1];
81 elseif strcmp(varargin{i},'classes')
82 classes = varargin{i+1};
83 eaten = [eaten, i, i+1];
84 end
85 end
86
87 global BuToolsCheckInput;
88
89 if isempty(BuToolsCheckInput)
90 BuToolsCheckInput = true;
91 end
92
93 if BuToolsCheckInput && ~CheckGenerator(Q)
94 error('FluidPrioQueue: Matrix Q is not a valid continuous time generator matrix!');
95 end
96
97 if BuToolsCheckInput && d<precision
98 error('FluidPrioQueue: The fluid service rate must be positive!');
99 end
100
101 if BuToolsCheckInput
102 for k=1:K
103 if ~all(all(R>=-precision))
104 error('FluidPrioQueue: The fluid arrival rate can not be negative!');
105 end
106 end
107 end
108
109 % Auxiliary functions
110 % ===================
111
112 % calculates the nth derivative of inv(v*R-Q), even if R contains zero elements
113 function dr = DReward (Q, R, n)
114 NQ = size(Q,1);
115 ix = (1:NQ);
116 ixz = ix(abs(diag(R))<=precision);
117 ixp = [ix(diag(R)>precision), ix(diag(R)<-precision)];
118 Nz = length(ixz);
119 Np = length(ixp);
120 Per = zeros(NQ);
121 for iw=1:Nz
122 Per(iw,ixz(iw))=1;
123 end
124 for iw=1:Np
125 Per(Nz+iw,ixp(iw))=1;
126 end
127 iPer = inv(Per);
128 Rp = R(ixp,ixp);
129 Qpp = Q(ixp,ixp);
130 Qpz = Q(ixp,ixz);
131 Qzp = Q(ixz,ixp);
132 Qzz = Q(ixz,ixz);
133 dXvn = (-1)^n * factorial(n) * inv(inv(Rp)*(-Qpp-Qpz*inv(-Qzz)*Qzp))^(n+1) * inv(Rp);
134 drpar = [inv(-Qzz)*Qzp*dXvn*Qpz*inv(-Qzz), inv(-Qzz)*Qzp*dXvn; dXvn*Qpz*inv(-Qzz), dXvn];
135 dr = iPer*drpar*Per;
136 end
137
138 function BPM = BusyPeriodRewardMoms (F, C, D, numOfMoms)
139 % block partitioning of the input
140 NF = size(F,1);
141 ix = (1:NF);
142 ixz = ix(abs(diag(C))<=precision);
143 ixp = ix(diag(C)>precision);
144 ixn = ix(diag(C)<-precision);
145 Nz = length(ixz);
146 Np = length(ixp);
147 Nn = length(ixn);
148 % permutation matrix that converts between the original and the partitioned state ordering
149 Per = zeros(NF);
150 for i=1:Nz
151 Per(i,ixz(i))=1;
152 end
153 for i=1:Np
154 Per(Nz+i,ixp(i))=1;
155 end
156 for i=1:Nn
157 Per(Nz+Np+i,ixn(i))=1;
158 end
159 iPer = inv(Per);
160 Fc = mat2cell(Per*F*iPer, [Nz, Np, Nn], [Nz, Np, Nn]);
161 [Fzz,Fpz,Fmz,Fzp,Fpp,Fmp,Fzm,Fpm,Fmm] = Fc{:};
162 Cm = C(ixn,ixn);
163 Cp = C(ixp,ixp);
164 Dm = D(ixn,ixn);
165 Dp = D(ixp,ixp);
166 Dz = D(ixz,ixz);
167
168 % detivatives of F(v)
169 Fppd = cell(1,numOfMoms+1);
170 Fpmd = cell(1,numOfMoms+1);
171 Fmpd = cell(1,numOfMoms+1);
172 Fmmd = cell(1,numOfMoms+1);
173 Fppd{1} = inv(Cp)*(Fpp+Fpz*inv(-Fzz)*Fzp);
174 Fpmd{1} = inv(Cp)*(Fpm+Fpz*inv(-Fzz)*Fzm);
175 Fmpd{1} = inv(-Cm)*(Fmp+Fmz*inv(-Fzz)*Fzp);
176 Fmmd{1} = inv(-Cm)*(Fmm+Fmz*inv(-Fzz)*Fzm);
177 for i=1:numOfMoms
178 dr = DReward(Fzz, Dz, i);
179 Fppd{i+1} = inv(Cp) * Fpz * dr * Fzp;
180 Fpmd{i+1} = inv(Cp) * Fpz * dr * Fzm;
181 Fmpd{i+1} = inv(-Cm) * Fmz * dr * Fzp;
182 Fmmd{i+1} = inv(-Cm) * Fmz * dr * Fzm;
183 if i==1
184 Fppd{i+1} = Fppd{i+1} - inv(Cp)*Dp;
185 Fmmd{i+1} = Fmmd{i+1} - inv(-Cm)*Dm;
186 end
187 end
188 Psi = FluidFundamentalMatrices(Fppd{1}, Fpmd{1}, Fmpd{1}, Fmmd{1}, 'P', precision);
189 BPM = cell(1,numOfMoms+1);
190 BPM{1} = Psi;
191 for i=1:numOfMoms
192 X = -Psi*Fmpd{i+1}*Psi + Fpmd{i+1};
193 for m=0:i-1
194 X = X + nchoosek(i,m) * ((Fppd{i-m+1} + Psi*Fmpd{i-m+1})*BPM{m+1} + BPM{m+1}*(Fmmd{i-m+1}+Fmpd{i-m+1}*Psi));
195 end
196 for l=1:i-1
197 for m=1:i-l
198 X = X + nchoosek(i,l)*nchoosek(i-l,m)*BPM{l+1}*Fmpd{i-l-m+1}*BPM{m+1};
199 end
200 end
201 BPM{i+1} = lyap(Fppd{1}+Psi*Fmpd{1}, Fmmd{1}+Fmpd{1}*Psi, X);
202 end
203 % re-order states
204 for i=1:length(BPM)
205 BPM{i} = iPer*[zeros(Nz,NF);zeros(Np,Nz+Np),BPM{i};zeros(Nn,NF)]*Per;
206 end
207 end
208
209 function [pr, Pn] = BusyPeriodRewardDistr (F, C, D, t)
210 % block partitioning of the input
211 NF = size(F,1);
212 ix = (1:NF);
213 ixz = ix(abs(diag(C))<=precision);
214 ixp = ix(diag(C)>precision);
215 ixn = ix(diag(C)<-precision);
216 Nz = length(ixz);
217 Np = length(ixp);
218 Nn = length(ixn);
219 % permutation matrix that converts between the original and the partitioned state ordering
220 Per = zeros(NF);
221 for i=1:Nz
222 Per(i,ixz(i))=1;
223 end
224 for i=1:Np
225 Per(Nz+i,ixp(i))=1;
226 end
227 for i=1:Nn
228 Per(Nz+Np+i,ixn(i))=1;
229 end
230 iPer = inv(Per);
231 Fc = mat2cell(Per*F*iPer, [Nz, Np, Nn], [Nz, Np, Nn]);
232 [Fzz,Fpz,Fmz,Fzp,Fpp,Fmp,Fzm,Fpm,Fmm] = Fc{:};
233 Cm = C(ixn,ixn);
234 Cp = C(ixp,ixp);
235 Dm = D(ixn,ixn);
236 Dp = D(ixp,ixp);
237 Dz = D(ixz,ixz);
238 % start erlangization
239 L = erlMaxOrder;
240 nu = L/t;
241 Z = inv(nu*Dz-Fzz);
242 Psie = FluidFundamentalMatrices (inv(Cp)*(Fpp-nu*Dp+Fpz*Z*Fzp), inv(Cp)*(Fpm+Fpz*Z*Fzm), inv(-Cm)*(Fmp+Fmz*Z*Fzp), inv(-Cm)*(Fmm-nu*Dm+Fmz*Z*Fzm), 'P', precision);
243 Pn = {Psie};
244 pr = Psie;
245 AM = inv(Cp)*(Fpp-nu*Dp+Fpz*Z*Fzp) + Psie*inv(-Cm)*(Fmp+Fmz*Z*Fzp);
246 BM = inv(-Cm)*(Fmm-nu*Dm+Fmz*Z*Fzm) + inv(-Cm)*(Fmp+Fmz*Z*Fzp)*Psie;
247 for n=1:L-1
248 CM = inv(Cp)*nu*Dp*Pn{n} + Pn{n}*inv(-Cm)*nu*Dm;
249 for i=1:n-1
250 CM = CM + Pn{i+1}*inv(-Cm)*Fmp*Pn{n-i+1};
251 end
252 CM = CM + inv(Cp)*Fpz*Z*(nu*Dz*Z)^n*Fzm - Psie*inv(-Cm)*Fmz*Z*(nu*Dz*Z)^n*Fzp*Psie;
253 if ~isempty(ixz)
254 for i=0:n-1
255 CM = CM + Pn{i+1}*inv(-Cm)*Fmz*Z*(nu*Dz*Z)^(n-i)*(Fzm+Fzp*Psie);
256 CM = CM + (inv(Cp)*Fpz+Psie*inv(-Cm)*Fmz)*Z*(nu*Dz*Z)^(n-i)*Fzp*Pn{i+1};
257 end
258 for i=1:n-1
259 for j=1:n-i
260 CM = CM + Pn{i+1}*inv(-Cm)*Fmz*Z*(nu*Dz*Z)^(n-i-j)*Fzp*Pn{j+1};
261 end
262 end
263 end
264 PM = lyap(AM, BM, CM);
265 Pn{n+1} = PM;
266 % accumulation
267 pr = pr + PM;
268 end
269 % re-order states
270 pr = iPer*[zeros(Nz,NF);zeros(Np,Nz+Np),pr;zeros(Nn,NF)]*Per;
271 for i=1:length(Pn)
272 Pn{i} = iPer*[zeros(Nz,NF);zeros(Np,Nz+Np),Pn{i};zeros(Nn,NF)]*Per;
273 end
274 % end of erlangization
275 end
276
277 % some preparation
278 pi = CTMCSolve(Q);
279 lambda = pi*R';
280 N = size(Q,1);
281
282 % step 2. calculate performance measures
283 % ======================================
284 Ret = {};
285 for k=classes
286 % step 1. solution of the workload process for fluid types having
287 % the same or higher priority
288 % ============================================================
289 [mass0, ini, Km, clo] = GeneralFluidSolve (Q, diag(sum(R(k:end,:),1))/d-eye(N), [], precision);
290 KN = size(Km,1);
291 clok = clo*diag(R(k,:)) / lambda(k);
292
293 % similarity transformation
294 Delta = diag(inv(-Km)*sum(clok,2));
295 K0 = inv(Delta)*Km*Delta;
296 K1 = inv(Delta)*clok;
297 kappa = ini*Delta;
298
299 Km = K0;
300 clok = K1;
301 ini = kappa;
302
303 if k<K
304 % step 4.3. calculate the performance measures
305 % ==========================================
306 argIx = 1;
307 while argIx<=length(varargin)
308 if any(ismember(eaten, argIx))
309 argIx = argIx + 1;
310 continue;
311 elseif strcmp(varargin{argIx},'stMoms')
312 % MOMENTS OF THE SOJOURN TIME
313 % ~~~~~~~~~~~~~~~~~~~~~~~~~~~
314 numOfSTMoms = varargin{argIx+1};
315 F = [Km, clok; zeros(N,KN), Q];
316 C = [eye(KN), zeros(KN,N); zeros(N,KN), diag(sum(R(k+1:end,:),1))/d-eye(N)];
317 D = [zeros(KN,KN+N); zeros(N,KN), eye(N)];
318 Tmp = BusyPeriodRewardMoms (F, C, D, numOfSTMoms);
319 inis = [ini, zeros(1,N)];
320 stMoms = zeros(1,numOfSTMoms);
321 for i=1:length(Tmp)-1
322 stMoms(i) = (-1)^i*sum(inis * Tmp{i+1});
323 end
324 Ret{end+1} = stMoms;
325 argIx = argIx + 1;
326 elseif strcmp(varargin{argIx},'flMoms')
327 % MOMENTS OF THE FLUID LEVEL
328 % ~~~~~~~~~~~~~~~~~~~~~~~~~~
329 % first compute moments right after fluid drop
330 % departures
331 numOfFLMoms = varargin{argIx+1};
332 F = [Km, clok; zeros(N,KN), Q];
333 C = [eye(KN), zeros(KN,N); zeros(N,KN), diag(sum(R(k+1:end,:),1))/d-eye(N)];
334 D = [zeros(KN,KN+N); zeros(N,KN), diag(R(k,:))];
335 FLDPn = BusyPeriodRewardMoms (F, C, D, numOfFLMoms);
336 for i=1:length(FLDPn)
337 X = FLDPn{i};
338 FLDPn{i} = X(:,KN+1:end);
339 end
340 inis = [ini, zeros(1,N)];
341 fldMoms = zeros(1,numOfFLMoms);
342 for i=1:length(FLDPn)
343 FLDPn{i} = inis*FLDPn{i};
344 if i==1
345 FLDPn{i} = FLDPn{i} + mass0*diag(R(k,:))/lambda(k);
346 end
347 fldMoms(i) = (-1)^(i-1)*sum(FLDPn{i});
348 end
349 % calculate moments in random point of time
350 FLPn = {pi};
351 flMoms = zeros(1,numOfFLMoms);
352 iTerm = inv(ones(N,1)*pi - Q);
353 for n=1:numOfFLMoms
354 sumP = sum(FLDPn{n+1}) + n*(-FLDPn{n} + FLPn{n}*diag(R(k,:))/lambda(k))*iTerm*R(k,:)';
355 P = sumP*pi + n*(-FLPn{n}*diag(R(k,:)) + FLDPn{n}*lambda(k))*iTerm;
356 FLPn{n+1} = P;
357 flMoms(n) = (-1)^n*sum(P);
358 end
359 Ret{end+1} = flMoms;
360 argIx = argIx + 1;
361 elseif strcmp(varargin{argIx},'stDistr')
362 % DISTRIBUTION OF THE SOJOURN TIME
363 % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
364 stCdfPoints = varargin{argIx+1};
365 F = [Km, clok; zeros(N,KN), Q];
366 C = [eye(KN), zeros(KN,N); zeros(N,KN), diag(sum(R(k+1:end,:),1))/d-eye(N)];
367 D = [zeros(KN,KN+N); zeros(N,KN), eye(N)];
368 inis = [ini, zeros(1,N)];
369 res = zeros(1,length(stCdfPoints));
370 for x=1:length(stCdfPoints)
371 Tmp = BusyPeriodRewardDistr (F, C, D, stCdfPoints(x));
372 res(x) = sum(mass0*diag(R(k,:))/lambda(k)) + sum(inis*Tmp);
373 end
374 Ret{end+1} = res;
375 argIx = argIx + 1;
376 elseif strcmp(varargin{argIx},'flDistr')
377 % DISTRIBUTION OF THE FLUID LEVEL
378 % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
379 flCdfPoints = varargin{argIx+1};
380 F = [Km, clok; zeros(N,KN), Q];
381 C = [eye(KN), zeros(KN,N); zeros(N,KN), diag(sum(R(k+1:end,:),1))/d-eye(N)];
382 D = [zeros(KN,KN+N); zeros(N,KN), diag(R(k,:))];
383 inis = [ini, zeros(1,N)];
384 res = zeros(1,length(flCdfPoints));
385 % resDep = zeros(1,length(flCdfPoints));
386 for x=1:length(flCdfPoints)
387 nu = erlMaxOrder/flCdfPoints(x);
388 [Tmp, Psix] = BusyPeriodRewardDistr (F, C, D, flCdfPoints(x));
389 Psiy = lambda(k)*nu*(mass0*diag(R(k,:))/lambda(k)+inis*Psix{1}(:,KN+1:end)) * inv(nu*diag(R(k,:))-Q);
390 for i=2:length(Psix)
391 Psiy = nu*(lambda(k)*inis*Psix{i}(:,KN+1:end) + Psiy*diag(R(k,:)))*inv(nu*diag(R(k,:))-Q);
392 end
393 res(x) = sum(Psiy);
394 % resDep(x) = sum(mass0*diag(R(k,:))/lambda(k)) + sum(inis*Tmp);
395 end
396 Ret{end+1} = res;
397 argIx = argIx + 1;
398 else
399 error (['FluidPrioQueue: Unknown parameter ' varargin{argIx}])
400 end
401 argIx = argIx + 1;
402 end
403 elseif k==K
404 % step 3. calculate the performance measures
405 % ==========================================
406 argIx = 1;
407 while argIx<=length(varargin)
408 if any(ismember(eaten, argIx))
409 argIx = argIx + 1;
410 continue;
411 elseif strcmp(varargin{argIx},'stMoms')
412 % MOMENTS OF THE SOJOURN TIME
413 % ~~~~~~~~~~~~~~~~~~~~~~~~~~~
414 numOfSTMoms = varargin{argIx+1};
415 stMoms = zeros(1,numOfSTMoms);
416 for i=1:numOfSTMoms
417 stMoms(i) = factorial(i) * sum(ini*inv(-Km)^(i+1)*clok);
418 end
419 Ret{end+1} = stMoms;
420 argIx = argIx + 1;
421 elseif strcmp(varargin{argIx},'stDistr')
422 % DISTRIBUTION OF THE SOJOURN TIME
423 % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
424 stCdfPoints = varargin{argIx+1};
425 res = zeros(1,length(stCdfPoints));
426 for x=1:length(stCdfPoints)
427 res(x) = sum(mass0*diag(R(k,:))/lambda(k)) + sum(ini*inv(-Km)*(eye(size(Km))-expm(Km*stCdfPoints(x)))*clok);
428 end
429 Ret{end+1} = res;
430 argIx = argIx + 1;
431
432 elseif strcmp(varargin{argIx},'flMoms')
433 % MOMENTS OF THE FLUID LEVEL
434 % ~~~~~~~~~~~~~~~~~~~~~~~~~~
435 % first compute moments right after fluid drop
436 % departures
437 numOfFLMoms = varargin{argIx+1};
438 Ret{end+1} = FluFluQueue(Q,diag(R(k,:)),0,d,false,'flMoms', numOfFLMoms);
439 argIx = argIx + 1;
440 elseif strcmp(varargin{argIx},'flDistr')
441 % DISTRIBUTION OF THE FLUID LEVEL
442 % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
443 % at departures:
444 flCdfPoints = varargin{argIx+1};
445 Ret{end+1} = FluFluQueue(Q,diag(R(k,:)),0,d,false,'flDistr', flCdfPoints);
446 argIx = argIx + 1;
447 else
448 error (['FluidPrioQueue: Unknown parameter ' varargin{argIx}])
449 end
450 argIx = argIx + 1;
451 end
452 end
453 end
454 varargout = Ret;
455end
the
Definition Station.m:245
Generated by doxygen 1.9.8