4 вариант 10

from sympy import solve
from sympy.abc import x
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

lam = 1/5
mu = 1/3
a = 1
b = 8
n = 10
def get_r(x):
    return  x + x**2 + x**3 + x**4 + x**5 + x**6 - lam / mu
r = fsolve(get_r, [0,1])[1]

def get_p_0(r, a):
  return (a/(1 - r) + (r**(a + 1) - r**(b + 1))/(1 - r)**2)**(-1)

def get_p_q():
  return get_p_0(r, a)

def get_p_n(r, n):
    return (get_p_0(r, a) * r**(n + 1))*((1 - r**b)/(1 - r))

def get_avg_queue(r, a, b):
    p_0 = get_p_0(r, a)
    temp1 = p_0/(1-r)
    temp3 = (r**2*(a*r**(a-1)*(1-r)-(1-r**a)))/(1-r)**2
    temp4 = p_0*(1-r**b)/(1-r) * (r**2 / (1-r)**2)
    avg_queue = temp1 * temp3 + temp4
    return avg_queue

def get_w(lam, r, a, b):
  z = get_avg_queue(r, a, b)
  return z/lam
def get_z0(r,a):
    return 1/(1-get_p_0(r, a))

def get_n(r,a):
    return 1/(get_z0(r,a)-1)

print ('r = ', r)
print('Cреднее время ожидания подъемника: ', get_w(lam, r, a, b))
print('Cреднее число человек, находящихся  в  подъемнике: ', get_n(r,a)-get_avg_queue(r, a, b))
print('Среднее число человек, ожидающих подъема: ', get_avg_queue(r, a, b))
print('Вероятность ожидания в очереди: ', 1 - get_p_0(r, a))
print('Доля времени, в течение которого подъемник используется: ', 1 - get_p_0(r, a))
print('Вероятность  того,  что  в  подъемнике  находится максимальное число человек: ',1 - get_p_0(r, a) - sum([get_p_n(r, n) for n in range (b)]))

lam_arr = []
p_q_arr = []
w_arr = []

def get_rb7(x):
    return  x + x**2 + x**3 + x**4 + x**5 + x**6 + x**7 - lam / mu
def get_rb8(x):
    return  x + x**2 + x**3 + x**4 + x**5 + x**6 + x**7 + x**8 - lam / mu
def get_rb9(x):
    return  x + x**2 + x**3 + x**4 + x**5 + x**6 + x**7 + x**8 + x**9 - lam / mu
def get_rb10(x):
    return  x + x**2 + x**3 + x**4 + x**5 + x**6 + x**7 + x**8 + x**9 + x**10 - lam / mu
def get_rb11(x):
    return  x + x**2 + x**3 + x**4 + x**5 + x**6 + x**7 + x**8 + x**9 + x**10 + x**11 - lam / mu
def get_rb12(x):
    return  x + x**2 + x**3 + x**4 + x**5 + x**6 + x**7 + x**8 + x**9 + x**10 + x**11 + x**12 - lam / mu

loc_r7 = fsolve(get_rb7, [0,1])[1]
loc_r8 = fsolve(get_rb8, [0,1])[1]
loc_r9 = fsolve(get_rb9, [0,1])[1]
loc_r10 = fsolve(get_rb10, [0,1])[1]
loc_r11 = fsolve(get_rb11, [0,1])[1]
loc_r12 = fsolve(get_rb12, [0,1])[1]

lam_arr.append(7)
lam_arr.append(8)
lam_arr.append(9)
lam_arr.append(10)
lam_arr.append(11)
lam_arr.append(12)

p_q_arr.append(1 - get_p_0(loc_r7, a) - sum([get_p_n(loc_r7, n) for n in range (7)]))
p_q_arr.append(1 - get_p_0(loc_r8, a) - sum([get_p_n(loc_r8, n) for n in range (8)]))
p_q_arr.append(1 - get_p_0(loc_r9, a) - sum([get_p_n(loc_r9, n) for n in range (9)]))
p_q_arr.append(1 - get_p_0(loc_r10, a) - sum([get_p_n(loc_r10, n) for n in range (10)]))
p_q_arr.append(1 - get_p_0(loc_r11, a) - sum([get_p_n(loc_r11, n) for n in range (11)]))
p_q_arr.append(1 - get_p_0(loc_r12, a) - sum([get_p_n(loc_r12, n) for n in range (12)]))

w_arr.append(get_w(lam, loc_r7, a, 7))
w_arr.append(get_w(lam, loc_r8, a, 8))
w_arr.append(get_w(lam, loc_r9, a, 9))
w_arr.append(get_w(lam, loc_r10, a, 10))
w_arr.append(get_w(lam, loc_r11, a, 11))
w_arr.append(get_w(lam, loc_r12, a, 12))

fig, ax = plt.subplots()
ax.set_xlabel('Максимальная вместимость подъемника')
ax.set_ylabel('Вероятность  того,  что  в  подъемнике  находится максимальное число человек')
ax.plot(lam_arr, p_q_arr, marker='o')
for i in range(len(p_q_arr)):
    ax.annotate(f'({lam_arr[i]}; {round(p_q_arr[i], 3)})', (lam_arr[i], p_q_arr[i]))
plt.show()

fig, ax = plt.subplots()
ax.set_xlabel('Максимальная вместимость подъемника')
ax.set_ylabel('Cреднее время ожидания подъемника')
ax.plot(lam_arr, w_arr, marker='o')
for i in range(len(p_q_arr)):
    ax.annotate(f'({lam_arr[i]}; {round(w_arr[i], 3)})', (lam_arr[i], w_arr[i]))
plt.show()