MAM (Matrix Analytic Methods)

The MAM solver analyzes Markovian open queueing systems using matrix-analytic methods. The core solver is built on top of the BuTools library [1], Q-MAM, SMCSolver, and MAMSolverresume [3]. Native algorithms for parametric decomposition of queueing networks are integrated with this solver.

MethodAlgorithmReference
defaultMatrix-analytic solution of structured QBDs[1]
exactExact matrix-analytic solution
dec.sourceDecomposition with source-like arrivals
dec.mmapDecomposition with Markovian arrival process (MMAP)
dec.poissonDecomposition with Poisson arrival flows
mnaMatrix Network Analyzer decomposition[2]
inapIterative Numerical Algorithm for product-forms
inapplusIterative Numerical Algorithm for product-forms (enhanced)

MAP/BMAP/1 Queue Solver

The MAM solver includes specialized support for MAP/BMAP/1 queues (Markovian Arrival Process with Batch Markovian Service). This solver models the queue as a GI/M/1 type Markov chain (skip-free to the right) and uses the ETAQA algorithm for efficient computation of stationary probabilities and queue length moments.

Key Characteristics:

  • MAP arrivals: Arrivals follow a Markovian Arrival Process, adding one customer at a time
  • BMAP service: Service completions can process batches of 1, 2, 3, ... customers
  • GI/M/1-type structure: Level increases by exactly 1 (arrivals), decreases by any amount (batch service)

Features:

  • Supports phase-type modulated arrival processes (MAP)
  • Handles batch service with arbitrary batch size distributions
  • Computes mean queue length, utilization, response time, and throughput
  • Returns aggregated stationary probabilities and R matrix

MATLAB Usage:

% Create 2-phase MAP for arrivals
C0 = [-3.0, 1.0; 0.5, -2.0];  % transitions without arrivals
C1 = [1.5, 0.5; 0.5, 1.0];    % arrival transitions

% Create BMAP for service with batch sizes 1 and 2
D0 = [-5.0, 1.0; 0.5, -4.5];  % transitions without service
D1 = [2.0, 0.5; 0.5, 2.0];    % batch size 1 service
D2 = [1.0, 0.5; 0.5, 1.0];    % batch size 2 service

% Solve MAP/BMAP/1
[QN, UN, RN, TN] = solver_mam_map_bmap_1({C0, C1}, {D0, D1, D2});

% Display results
fprintf('Mean queue length: %.4f\n', QN);
fprintf('Utilization: %.4f\n', UN);
fprintf('Mean response time: %.4f\n', RN);
fprintf('Throughput: %.4f\n', TN);

Python Usage:

import numpy as np
from line_solver.api_native.mam import solver_mam_map_bmap_1

# Create 2-phase MAP for arrivals
C0 = np.array([[-3.0, 1.0], [0.5, -2.0]])
C1 = np.array([[1.5, 0.5], [0.5, 1.0]])

# Create BMAP for service with batch sizes 1 and 2
D0 = np.array([[-5.0, 1.0], [0.5, -4.5]])
D1 = np.array([[2.0, 0.5], [0.5, 2.0]])
D2 = np.array([[1.0, 0.5], [0.5, 1.0]])

# Solve MAP/BMAP/1
result = solver_mam_map_bmap_1(C0, C1, [D0, D1, D2])

print(f"Mean queue length: {result.mean_queue_length:.4f}")
print(f"Utilization: {result.utilization:.4f}")
print(f"Mean response time: {result.mean_response_time:.4f}")
print(f"Throughput: {result.throughput:.4f}")

BMAP/MAP/1 Queue Solver

The MAM solver includes specialized support for BMAP/MAP/1 queues (Batch Markovian Arrivals with Markovian Service). This solver models the queue as an M/G/1 type Markov chain (skip-free to the left) and uses the ETAQA algorithm for efficient computation of stationary probabilities and queue length moments.

Key Characteristics:

  • BMAP arrivals: Batch arrivals with sizes 1, 2, 3, ... following a Batch Markovian Arrival Process
  • MAP service: Service completions follow a Markovian Arrival Process, serving one customer at a time
  • M/G/1-type structure: Level increases by any amount (batch arrivals), decreases by exactly 1 (single service)

Features:

  • Supports batch arrivals with arbitrary batch size distributions
  • Handles phase-type modulated service processes (MAP)
  • Computes mean queue length, utilization, response time, and throughput
  • Returns aggregated stationary probabilities and G matrix

MATLAB Usage:

% Create BMAP for arrivals with batch sizes 1 and 2
D0 = [-2.8, 1.6; 1.6, -2.8];  % transitions without arrivals
D1 = [0.5, 0.2; 0.2, 0.5];    % batch size 1 arrivals
D2 = [0.3, 0.2; 0.2, 0.3];    % batch size 2 arrivals

% Create 2-phase MAP for service
S0 = [-5.0, 1.0; 0.5, -4.5];  % transitions without service
S1 = [2.5, 1.5; 1.5, 2.5];    % service transitions

% Solve BMAP/MAP/1
[QN, UN, RN, TN] = solver_mam_bmap_map_1({D0, D1, D2}, {S0, S1});

% Display results
fprintf('Mean queue length: %.4f\n', QN);
fprintf('Utilization: %.4f\n', UN);
fprintf('Mean response time: %.4f\n', RN);
fprintf('Throughput: %.4f\n', TN);

Python Usage:

import numpy as np
from line_solver.api_native.mam import solver_mam_bmap_map_1

# Create BMAP for arrivals with batch sizes 1 and 2
D0 = np.array([[-2.8, 1.6], [1.6, -2.8]])
D1 = np.array([[0.5, 0.2], [0.2, 0.5]])
D2 = np.array([[0.3, 0.2], [0.2, 0.3]])

# Create 2-phase MAP for service
S0 = np.array([[-5.0, 1.0], [0.5, -4.5]])
S1 = np.array([[2.5, 1.5], [1.5, 2.5]])

# Solve BMAP/MAP/1
result = solver_mam_bmap_map_1([D0, D1, D2], S0, S1)

print(f"Mean queue length: {result.mean_queue_length:.4f}")
print(f"Utilization: {result.utilization:.4f}")
print(f"Mean response time: {result.mean_response_time:.4f}")
print(f"Throughput: {result.throughput:.4f}")

Example

This example demonstrates solving a simple M/M/1 queue using the MAM solver. The unique feature shown here is the iterative matrix-analytic method which converges to the solution (note the "Iterations: 4" in the output).

% Create a simple M/M/1 queue
model = Network('M/M/1 Example');
source = Source(model, 'Source');
queue = Queue(model, 'Queue1', SchedStrategy.FCFS);
sink = Sink(model, 'Sink');

jobclass = OpenClass(model, 'Class1');
source.setArrival(jobclass, Exp(0.9));
queue.setService(jobclass, Exp(1.0));

model.link(Network.serialRouting(source, queue, sink));

% Solve with MAM (Matrix Analytic Methods)
solver = MAM(model);

% Display average performance metrics
MAM(model).avgTable()

Output:

Default method: using dec.source
MAM analysis [method: default/dec.source, lang: matlab] completed in 0.123s. Iterations: 4.

ans =
  2×8 table
    Station    JobClass    QLen    Util    RespT    ResidT    ArvR    Tput
    _______    ________    ____    ____    _____    ______    ____    ____
    Source      Class1      0        0       0         0        0     0.9
    Queue1      Class1      9      0.9      10        10      0.9     0.9 
import jline.lang.*;
import jline.lang.constant.SchedStrategy;
import jline.lang.nodes.*;
import jline.lang.processes.Exp;
import jline.solvers.NetworkAvgTable;
import jline.solvers.mam.MAM;

public class MAMExample {
    public static void main(String[] args) {
        // Create a simple M/M/1 queue
        Network model = new Network("M/M/1 Example");

        Source source = new Source(model, "Source");
        Queue queue = new Queue(model, "Queue1", SchedStrategy.FCFS);
        Sink sink = new Sink(model, "Sink");

        OpenClass jobclass = new OpenClass(model, "Class1");
        source.setArrival(jobclass, Exp.fitRate(0.9));
        queue.setService(jobclass, Exp.fitRate(1.0));

        model.link(Network.serialRouting(source, queue, sink));

        // Solve with MAM (Matrix Analytic Methods)
        MAM solver = new MAM(model);

        // Display average performance metrics
        NetworkAvgTable avgTable = solver.avgTable;
        System.out.println(avgTable);
    }
}

Output:

Default method: using dec.source
MAM analysis [method: default/dec.source, lang: java] completed. Iterations: 4.

Station    JobClass    QLen    Util    RespT    ResidT    ArvR    Tput
Source     Class1      0.0     0.0     0.0      0.0       0.0     0.9
Queue1     Class1      9.0     0.9     10.0     10.0      0.9     0.9
from line_solver import *

# Create a simple M/M/1 queue
model = Network('M/M/1 Example')

source = Source(model, 'Source')
queue = Queue(model, 'Queue1', SchedStrategy.FCFS)
sink = Sink(model, 'Sink')

jobclass = OpenClass(model, 'Class1')
source.setArrival(jobclass, Exp(0.9))
queue.setService(jobclass, Exp(1.0))

model.link(Network.serialRouting(source, queue, sink))

# Solve with MAM (Matrix Analytic Methods)
solver = MAM(model)

# Display average performance metrics
print(solver.avg_table)

Output:

Default method: using dec.source
MAM analysis [method: default/dec.source, lang: python] completed. Iterations: 4.

  Station  JobClass  QLen  Util  RespT  ResidT  ArvR  Tput
0  Source    Class1   0.0   0.0    0.0     0.0   0.0   0.9
1  Queue1    Class1   9.0   0.9   10.0    10.0   0.9   0.9

References

  1. Horváth, G., & Telek, M. (2017). BuTools: Burstiness and heavy traffic modeling. IEEE Communications Surveys & Tutorials.
  2. Li, Z., & Casale, G. (2024). Matrix Network Analyzer: A new decomposition algorithm for phase-type queueing networks. Proc. of ACM/SPEC ICPE, 77-88.
  3. Riska, A., & Smirni, E. (2003). ETAQA: An Efficient Technique for the Analysis of QBD-Processes by Aggregation. Performance Evaluation, 54(2):151-177. MAMSolver