← Back to Tutorials Prev: Tutorial 09 Next: Tutorial 11

Tutorial 10: Departure Process Analysis

This example illustrates LINE's support for extracting simulation data about particular events, such as departures from a queue. We compute the squared coefficient of variation of inter-departure times and compare with Marshall's exact formula.

M/E/1 Queue Network Diagram 
% Example 9: Studying a departure process
[model,source,queue,sink,oclass] = gallery_merl1;
%% Block 4: solution
solver = CTMC(model,'cutoff',150,'seed',23000);

sa = solver.sampleSysAggr(5e3);
ind = model.getNodeIndex(queue);

filtEvent = cellfun(@(c) c.node == ind && ...
    (isequal(c.event, EventType.DEP) || (isnumeric(c.event) && c.event == 2)), sa.event);
interDepTimes = diff(cellfun(@(c) c.t, {sa.event{filtEvent}}));

% Estimated squared coeff. of variation of departures
SCVdEst = var(interDepTimes)/mean(interDepTimes)^2

util = solver.getAvgUtil();
util = util(queue);
avgWaitTime = solver.getAvgWaitT();
avgWaitTime = avgWaitTime(queue);
SCVa = source.getArrivalProcess(oclass).getSCV();
svcRate = queue.getServiceProcess(oclass).getRate();
SCVs = queue.getServiceProcess(oclass).getSCV();

% Marshall's exact formula
SCVd = SCVa + 2*util^2*SCVs - 2*util*(1-util)*svcRate*avgWaitTime

fprintf('Simulated SCV:   %.6f\n', SCVdEst);
fprintf('Theoretical SCV: %.6f\n', SCVd);
from line_solver import *
import numpy as np

model, source, queue, sink, oclass = gallery_merl1()

# Block 4: solution
solver = CTMC(model, cutoff=150, seed=23000)
sa = solver.sample_sys_aggr(5000)
ind = model.get_node_index(queue)

# Filter departure events
dep_times = []
for event in sa.event:
    if event.node == ind and event.event == "DEP":
        dep_times.append(event.t)

if len(dep_times) > 1:
    inter_dep_times = np.diff(dep_times)
    # Estimated SCV of departures
    scv_d_est = np.var(inter_dep_times) / np.mean(inter_dep_times)**2
    print(f"Simulated SCV of departures: {scv_d_est:.6f}")
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import jline.examples.java.GettingStarted;
import jline.io.Ret.SampleResult;
import jline.lang.constant.EventType;
import jline.solvers.ctmc.CTMC;
import jline.util.matrix.Matrix;

Object[] components = GettingStarted.gallery_merl1();
Network model = (Network) components[0];
Source source = (Source) components[1];
Queue queue = (Queue) components[2];
OpenClass jobclass = (OpenClass) components[4];

CTMC solver = new CTMC(model, "cutoff", 150, "seed", 23000);
SampleResult sa = solver.sampleSys(5000);
int queueIndex = model.getNodeIndex(queue);

Matrix events = sa.event;
Matrix times = sa.t;
List<Double> departureTimes = new ArrayList<>();
for (int i = 0; i < events.getNumRows(); i++) {
    int eventType = (int) events.get(i, 0);
    int nodeIndex = (int) events.get(i, 1);
    if (eventType == EventType.DEP.ordinal() && nodeIndex == queueIndex)
        departureTimes.add(times.get(i, 0));
}

double[] interDepartureTimes = new double[departureTimes.size() - 1];
for (int i = 1; i < departureTimes.size(); i++)
    interDepartureTimes[i - 1] = departureTimes.get(i) - departureTimes.get(i - 1);
double mean = Arrays.stream(interDepartureTimes).average().orElse(0.0);
double variance = Arrays.stream(interDepartureTimes)
        .map(x -> Math.pow(x - mean, 2)).average().orElse(0.0);
double scvEstimate = variance / (mean * mean);
System.out.println("Simulated SCV: " + scvEstimate);
#include "line/solvers/ctmc/solver_ctmc_sample.h"
#include "line/solvers/ctmc/solver_ctmc_getters.h"

Network model("M/E/1");
Source source(model, "Source");
Queue queue(model, "Queue", SchedStrategy::FCFS);
Sink sink(model, "Sink");
OpenClass jobclass(model, "Class1");
source.set_arrival(jobclass, Exp(1.0));
queue.set_service(jobclass, Erlang(4.0, 2)); // mean 0.5
Routing P;
P.set(source, queue, 1.0);
P.set(queue, sink, 1.0);
model.link(P);

ctmc::CtmcOptions opt;
opt.cutoff = 150;
const auto sa = ctmc::solver_ctmc_sample_sys(
        model.get_struct(), opt, 5000, 23000);
const auto gen = ctmc::ctmc_get_infgen(model.get_struct(), sa.chain);
const std::size_t queueNode = model.get_struct().station_to_node[1];

std::vector<double> departureTimes;
for (std::size_t i = 0; i < sa.event.size(); ++i) {
    const std::size_t event = sa.event[i];
    if (event == static_cast<std::size_t>(-1) || event >= gen.sync.size()) continue;
    const auto& active = gen.sync[event].active;
    if (active.node == queueNode && active.event == lang::EventType::DEP)
        departureTimes.push_back(sa.t[i]);
}
std::vector<double> interDepartureTimes;
for (std::size_t i = 1; i < departureTimes.size(); ++i)
    interDepartureTimes.push_back(departureTimes[i] - departureTimes[i - 1]);
double mean = 0.0, secondMoment = 0.0;
for (double x : interDepartureTimes) { mean += x; secondMoment += x * x; }
mean /= interDepartureTimes.size();
secondMoment /= interDepartureTimes.size();
const double scvEstimate = (secondMoment - mean * mean) / (mean * mean);

Expected Output

=== Departure Process Analysis Results ===
Simulated SCV of departures:   0.899211
Theoretical SCV (Marshall):    0.875000
Relative error:                2.77%
← Back to Tutorials Prev: Tutorial 09 Next: Tutorial 11