Лаба 4 / 21_var

import sys
import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import erlang, expon

warnings.filterwarnings("ignore")


def plot_erlang_distribution(k, lambda_, figure):
    # Определение границ графика:
    temp = erlang.rvs(a=k, scale=1/lambda_, size=10000)
    x_min, x_max = min(temp), max(temp)

    # Построение графика
    x = np.linspace(x_min, x_max, 1000)
    y = erlang.pdf(x, k, scale=1/lambda_)
    plt.figure(figure)
    plt.plot(x, y, label=fr'$k = {k}$' + "\n" + fr'$\lambda = {lambda_}$', color='black')
    plt.grid()
    plt.legend()
    plt.title(f'Плотность эрланговского распределения')

def plot_expon_distribution(lambda_, figure):
    # Определение границ графика:
    temp = expon.rvs(scale=1/lambda_, size=10000)
    x_min, x_max = min(temp), max(temp)

    # Построение графика
    x = np.linspace(x_min, x_max, 1000)
    y = expon.pdf(x, scale=1/lambda_)
    plt.figure(figure)
    plt.plot(x, y, label=fr'$\lambda = {lambda_}$', color='black')
    plt.grid()
    plt.legend()
    plt.title(f'Плотность экспоненциального распределения распределения')

def model(mu, k_, l, is_test):
    # Основные параметры
    t_income = (
        lambda: np.random.exponential(1 / l)
        if is_test
        else erlang.rvs(a=k_, scale=1 / (l * k_))
    )
    t_outcome = lambda: np.random.exponential(scale=1 / mu)
    ts = 0
    busy = False
    t_in = t_income()
    t_ins = []
    t_out = t_in
    t_outs = []
    n, k, m = 0, 0, 0
    t_mean_old = sys.float_info.max
    counter = 0

    # Процесс моделирования
    while True:
        counter += 1
        # Проверяем наступило ли время освобождения ОУ
        if t_in <= t_out:
            # Записываем в системное время момент поступления заявки
            ts = t_in
            n += 1  # Увеличиваем счетчик поступивших заявок
            # Сохраняем в список момент поступления заявки
            t_ins.append(t_in)
            # Система занята?
            if busy:
                m += 1
            else:
                busy = True
                t_out = ts + t_outcome()
            # Вычисляем следующий момент поступления заявки
            t_in = ts + t_income()
        else:
            # Записываем в системное время момент освобождения ОУ
            ts = t_out
            k += 1  # Увеличиваем счетчик обслуженных заявок
            # Сохраняем в список момент освобождения ОУ
            t_outs.append(t_out)
            # В буфере больше 0 заявок?
            if m > 0:
                m -= 1
                t_out = ts + t_outcome()
            else:
                busy = False
                t_out = t_in
        # Рассматриваем возможность выхода каждые 1000 итераций
        if counter % 1000 == 0:
            t_mean_new = np.mean([t_outs[i] - t_ins[i] for i in range(len(t_outs))])
            if abs((t_mean_new - t_mean_old) / t_mean_old) < 0.0001:
                return t_mean_new
            t_mean_old = t_mean_new

def testing_model(mu, k, is_test, figure):
    # Параметры для тестирования модели
    lambdas = np.linspace(0.1, 1, 10) * mu
    T_theoretical = []
    T_experimental = []

    # Сбор практических данных
    for l in lambdas:
        print(f'{l = :.1f}, {mu = }')
        T_experimental.append(model(mu, k, l, is_test))

    # Сбор тестовых данных + вывод их на график
    if is_test:
        for l in lambdas:
            ro = l / mu
            N_mean = ro / (1 - ro)
            T_theoretical.append(N_mean / l)

        plt.figure(figure)
        plt.plot(lambdas[:9], T_theoretical[:9], color="blue", label="Теоретическое")
        plt.plot(lambdas[:9], T_experimental[:9], color="orange", label="Экспериментальное")
        plt.legend()
        plt.title("Оценка результатов моделирования")
        plt.ylabel("Ср. время в системе")
        plt.xlabel("Интенсивность входного потока $(\lambda)$")
        plt.grid()

    return T_experimental

if __name__ == "__main__":
    # Параметры, заданные по варианту:
    mu = 1
    k = 2

    # Построение графиков функций распределений
    plot_erlang_distribution(k=k, lambda_=mu, figure=1)
    plot_expon_distribution(lambda_=mu, figure=2)

    # Запуск моделирования на тестовом случае
    print('Моделирование тестового случая: ')
    testing_model(mu, k, is_test=True, figure=3)

    # plt.show() # Импровизированная точка останова

    # Запуск моделирования на основном случае
    print('\nМоделирование основного случая:')
    index = [rf"{i:.1f}mu" for i in np.linspace(0.1, 1, 10)]    # Индексы для создания df
    arr = testing_model(mu, k, is_test=False, figure=None)
    df = pd.DataFrame(arr, index=index, columns=["T"])

    print(df.T)

    # Построение графика экспериментальных данных
    plt.figure(4)
    plt.plot(df[:9], color='black')
    plt.grid()
    plt.title("Моделирование для экспериментальной зависимости")
    plt.ylabel("Ср. время в системе")
    plt.xlabel("Интенсивность входного потока $(\lambda)$")

    plt.show()