qbd_raprap1.m

%{ @file qbd_raprap1.m
 %  @brief Solves a RAP/RAP/1 queue using QBD methods
 %
 %  @author LINE Development Team
%}
 
%{
 % @brief Analyzes RAP/RAP/1 queue using Quasi-Birth-Death process
 %
 % @details
 % This function solves a RAP/RAP/1 queue (Rational Arrival Process) using
 % QBD methods, computing throughput, queue length, utilization, and other
 % performance metrics.
 %
 % The two RAPs are INDEPENDENT of each other, so the QBD phase space is the
 % product of the two phase spaces and the blocks factor as Kronecker
 % products. This function is the thin product-space wrapper; the analysis
 % itself is the block-level core qbd_rap.m, which solves an arbitrary QBD
 % with RAP components. A model whose arrival process and sequence of
 % service times share a phase space, and are therefore cross-correlated,
 % has no such product structure and must call qbd_rap directly.
 %
 % @par References:
 % N. G. Bean and B. F. Nielsen, "Quasi-Birth-and-Death Processes with
 % Rational Arrival Process Components", Stochastic Models, 26(3), 2010,
 % pp. 309-334. The analysis rests on the prediction-process interpretation
 % of a RAP due to Asmussen and Bladt, which is what allows a QBD argument
 % to be carried over to matrices that are not nonnegative. The same
 % prediction process underlies the conditional-vector RAP sampler in
 % rap_sample.m.
 %
 % @par Phase ordering:
 % The QBD phase is the pair (arrival phase, service phase) laid out with
 % the ARRIVAL phase major and the service phase minor, i.e. the phase index
 % is (a-1)*ns + s. That is the ordering produced by kron(RAPa, eye(ns)) and
 % kron(eye(na), RAPs), and pqueue is indexed by it downstream in
 % solver_mam_basic, so the two Kronecker factors must not be swapped.
 %
 % @par Syntax:
 % @code
 % [XN, QN, UN, pqueue, R, eta, G, B, L, F] = qbd_raprap1(RAPa, RAPs)
 % [XN, QN, UN, pqueue, R, eta, G, B, L, F] = qbd_raprap1(RAPa, RAPs, util)
 % @endcode
 %
 % @par Parameters:
 % <table>
 % <tr><th>Name<th>Description
 % <tr><td>RAPa<td>Arrival process (RAP)
 % <tr><td>RAPs<td>Service process (RAP)
 % <tr><td>util<td>(Optional) Target utilization to scale service rate
 % </table>
 %
 % @par Returns:
 % <table>
 % <tr><th>Name<th>Description
 % <tr><td>XN<td>System throughput
 % <tr><td>QN<td>Mean queue length
 % <tr><td>UN<td>Utilization
 % <tr><td>pqueue<td>Queue length distribution
 % <tr><td>R<td>Rate matrix R
 % <tr><td>eta<td>Caudal characteristic
 % <tr><td>G<td>Rate matrix G
 % <tr><td>B<td>Backward transition block
 % <tr><td>L<td>Local transition block
 % <tr><td>F<td>Forward transition block
 % </table>
%}
function [XN,QN,UN,pqueue,R,eta,G,B,L,F]=qbd_raprap1(RAPa,RAPs,util)
% [XN,QN,UN,PQUEUE,R,ETA]=QBD_RAPRAP1(RAPA,RAPS,UTIL)
 
na = length(RAPa{1});
ns = length(RAPs{1});
 
if nargin>=3 %exist('util','var')
    RAPs = map_scale(RAPs,util/map_lambda(RAPa));
end
 
% QBD blocks of the RAP/RAP/1 queue: a level is the number in system, a
% phase is the (arrival,service) RAP phase pair, arrival phase major. The
% level rises on an arrival (F), falls on a service completion (B), and the
% two RAPs evolve independently between level changes (L, the Kronecker sum
% of the two hidden generators). At level 0 the queue is empty, so no
% service completion can occur and only the arrival RAP evolves: the
% boundary local block is B1 = kron(Ca,I) and the boundary up-block is F.
F = kron(RAPa{2},eye(ns));           % arrivals, level up
L = kron(RAPa{1},eye(ns)) + kron(eye(na),RAPs{1});
B = kron(eye(na),RAPs{2});           % service completions, level down
B1 = kron(RAPa{1},eye(ns));
 
% Theorem 7 of Bean and Nielsen (2010), in qbd_rap.m: G from the quadratic
% matrix equation, U = L + F*G, R = F*inv(-U), and pi0 from the boundary
% equation pi0*(B1 + R*B) = 0 normalised by pi0*inv(I-R)*e = 1.
[~,~,R,G,~,eta,~,pi0] = qbd_rap(F, L, B, F, B1, 0);
 
% Level series. The truncation rule is the one of QBD_pi with MaxNumComp
% 100: accumulate levels until the mass reaches 1-1e-10, capped at 101
% level vectors. qbd_rap returns the exact mean queue length in closed
% form, but this function keeps reporting the truncated series because
% pqueue is the documented return value and downstream consumers (and the
% JAR and Python ports) must cut the tail at exactly the same point,
% otherwise the means disagree at the 1e-9 level.
maxNumComp = 100;
pqueue = pi0;
sumpi = sum(pi0);
numit = 1;
while sumpi < 1-1e-10 && numit < 1+maxNumComp
    pqueue(numit+1,1:(na*ns)) = pqueue(numit,:)*R;
    numit = numit+1;
    sumpi = sumpi + sum(pqueue(numit,:));
end
 
numLevels = size(pqueue,1);
levelProb = sum(pqueue,2);
QN = (0:(numLevels-1))*levelProb;
 
if na == 1 && ns == 1
    UN = 1 - pqueue(1);
else
    UN= 1 - sum(pqueue(1,:));
end
XN=map_lambda(RAPa);
end