Gerasimov's residue (closed-form) normalizing constant, generalized to R classes. More...

#include <algorithm>
#include <cmath>
#include <limits>
#include <type_traits>
#include <vector>
#include "line/api/pfqn/pfqn_asympt_common.h"
#include "line/api/pfqn/pfqn_ca.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for pfqn_gerasimov.h:

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
NcResult< T >	line::pfqn::pfqn_gerasimov (const Matrix< T > &L, const std::vector< int > &N, const Matrix< T > &Z, double tol=1e-12, std::size_t maxterms=200000)
	Exact normalizing constant of a closed multiclass product-form network by ITERATED RESIDUES of its rational generating function, one class at a time.
template<class T>
NcResult< T >	line::pfqn::pfqn_gerasimov (const Matrix< T > &L, const std::vector< int > &N)
	Delay-free overload.

Detailed Description

Gerasimov's residue (closed-form) normalizing constant, generalized to R classes.

Templated port of matlab/src/api/pfqn/pfqn_gerasimov.m and jar/src/main/java/jline/api/pfqn/nc/Pfqn_gerasimov.java.

A. I. Gerasimov, "On Normalizing Constants in Multiclass Queueing Networks", Operations Research 43(4):704-711, 1995, evaluates

G(N_1,...,N_R) = (2 pi i)^-R int_G1 ... int_GR prod_s z_s^{N_s-1} prod_i (1 - sum_s x_is/z_s)^-1

by residues, and gives the resulting CLOSED FORM only for R = 1 (Thm 1-2) and R = 2 (Thm 3 for simple poles, Thm 4 for multiple ones), stating that "for three or more classes of customers, the normalizing constants can be found by numerical methods". This header implements the residue elimination itself, so the closed form is produced for ANY R; at R = 2 it reproduces Thm 3/4 term by term.

Written as a coefficient of the u_s = 1/z_s series,

G(N) = [prod_s u_s^{N_s}] exp(sum_s Z_s u_s) prod_i (1 - sum_s x_is u_s)^-1,

every factor is AFFINE in u, so singling out u_r gives f = A - B u_r with A affine in the surviving variables. Partial fractions in u_r,

[u_r^n] prod_j (A_j - B_j u_r)^-m_j = sum_j sum_{k=0}^{m_j-1} (-B_j)^-k C(n+m_j-k-1,n) B_j^n A_j^-(n+m_j-k) [t^k] prod_{l!=j} (C_jl - B_l t)^-m_l, C_jl = (A_l B_j - B_l A_j)/B_j,

map a sum of products of affine powers into another one with one variable fewer, and R-1 such steps leave a univariate coefficient extraction. At R = 2 the single step returns one term per station i, with outer factor x_i2^{N_2+M-1}/prod_{k!=i}(x_i2-x_k2), a pole of order N_2+1 at x_i1 and simple poles at the paper's z_1ik = (x_k1 x_i2 - x_i1 x_k2)/(x_i2 - x_k2): exactly Thm 3, with the multiple poles of Thm 4 (his xi_i < M) handled by the same step. Tied x_i2, vanishing x_i2 and identical station rows, all outside the paper's hypotheses, are ordinary cases here.

Cost. Let M be the number of stations and order the populations N_(1) <= ... <= N_(R). The first elimination turns the single input term into M, and every later one multiplies the count by C(S+M-1,M-1) + M-1, where S is the total population already eliminated: a pole of order S+1 has to be differentiated against the M-1 remaining ones. The innermost extraction then convolves M series of length N_(1). Hence R = 1 costs O(M N), Buzen's own cost; R = 2 costs O(M^2 N_(1)^2), INDEPENDENT OF N_(2); and R >= 3 costs the same times prod_{r=3}^{R} C(N_(r)+M-1, M-1). The R = 2 line is the reason to reach for this method: a population removed by residues enters only as a pole ORDER, i.e. through binomial coefficients, so it costs nothing at all. For R >= 3 the term count is polynomial in the populations of degree (M-1)(R-2) and exponential in R, which is why the paper stops at two classes and why maxterms exists.

Conditioning. The class left for the innermost extraction is the one with the SMALLEST population, because that population is the degree the final, sign-indefinite series is carried to, whereas the eliminated ones cancel nothing; basing on N = 100 instead of N = 6 in one four-station model cost 39 nats of lG in double.

Arithmetic: RATIONAL. Every operation is a product, a quotient or a binomial coefficient, so the exact backend runs the residue expansion with NO cancellation at all. That is the reason to have this method: in double the alternating sum loses accuracy (~1e-15 typically, ~1e-9 on ill-conditioned demand matrices), and the same expansion in exact arithmetic returns G to the last digit. Two affine forms count as one pole when they are proportional, which in the exact backend means proportional exactly, and in a floating one within the relative tolerance.

Definition in file pfqn_gerasimov.h.

Namespaces

Functions

Detailed Description