Объектная среда для вероятностного анализа в системе MATLAB

Вариант семинара (практического занятия) по вероятностному анализу систем в "Объектной среде для вероятностного анализа в системе MATLAB" с использованием электронных формул.

Методическая часть

Технические сложности вероятностного анализа систем связаны, в первую очередь, с общей
методической проблемой инженерной науки. Закономерности функционирования объектов проектирования (исследования) выражаются детерминированными соотношениями и алгоритмами, а влияние случайных факторов исследуется в рамках специально создаваемой вероятностной модели, ориентированной на математический аппарат теории вероятностей. Компьютеризация инженерных расчетов значительно продвинула детерминированный анализ, тогда как аппарат прикладной теории вероятностей остается прежним, заложенным Я. Бернулли, изложенным в учебниках Вентцель Е.С. и др. Адекватная компьютеризация, отражающая истинную, то есть объектную, а не алгоритмическую природу моделируемых процессов, открывает возможность инкапсуляции в структуры объектов их случайной природы (законов распределения). Переопределение операций по правилам теории вероятностей и выполнение расчетов по моделям, отражающим как физическую, так и случайную природу реальных процессов, даст представление о возможном поведении объекта исследования, которое не нуждается в дополнительном (искусственном) вероятностном анализе.

 "Объектная среда для вероятностного анализа в системе MATLAB" базируется на понятии "электронная формула". Электронные формулы в рабочей области MATLAB синтаксически оформляются как вызовы функций, имеющих содержательное толкование в контексте теории вероятностей и выполняющие действия, зависящие от комбинации аргументов. Аргументы электронных формул, как правило, различаются по смыслу, будучи объектами определенных классов: событий, случайных величин (СВ), геометрических фигур и т.д. Вычисления при исполнении функций выполняют объекты во взаимодействии между собой. Математический аппарат теории вероятностей реализован иерархией классов, исходные данные вероятностных задач принимаются конструкторами классов, создающих объекты в контексте задачи, электронные формулы управляют взаимодействием объектов и формируют результат, соответствующий комбинации аргументов.

Объектно-ориентированная технология вероятностных расчетов разработана для решения разнообразных задач исследования операций, связанных с анализом случайных событий, алгеброй случайных величин и функциональными отношениями между ними, системами случайных векторов. Интерпретирующие возможности объектной среды создаются целесообразной организацией системы классов и функций, построенной на теоретических основах прикладной теории вероятностей.

Включает в себя следующие разделы:

 -      вероятности сложных событий;

-      формула Байеса;

-      формула Бернулли;

-      дискретные случайные величины;

-      базовый класс непрерывных случайных величин;

-      показательный закон и связанные с ним классы распределений;

-      нормальный закон распределения;

-      класс двумерных нормально распределенных векторов;

-      функции случайных величин;

-      мультипликативные функции случайных величин;

-      аддитивные функции случайных величин.

Системные требования и установка

MATLAB версии не ниже R2019a. Для установки из Интернета (рекомендуемый способ) выполните в MATLAB команду eval(webread('https://exponenta.ru/install/m4')). Для оффлайн-установки скачайте и запустите файл ProbabilityAnalysisM4.mltbx (https://github.com/ETMC-Exponenta/ProbabilityAnalysis-M4/raw/master/ProbabilityAnalysisM4.mltbx). Для обновления библиотеки выполните в MATLAB команду m4update.

Пример применения электронной формулы

Если X – случайная величина с известным законом распределения, G – геометрическая фигура или массив непересекающихся фигур, электронная формула [P, p] = Ver(X, G) вычисляет вероятность попадания P случайной величины X в область G. Объект X создается конструктором одного из классов случайных величин с соответствующими параметрами, например, X = Norm_2([1; 2/3], [3, 2.4], 0.6) – это объект нормального распределения на плоскости с математическими ожиданиями (МО) по координатам, равными 1 и 2/3, среднеквадратическими отклонениями (СКО) 3 и 2,4, коэффициентом корреляции между координатами 0,6. Объект G может быть создан конструктором кругов заданного радиуса с центром в определенной точке G = Circ(1.75, [0; 1.75]), прямоугольников G = Rect(5, 3.5) или других фигур, может быть также результатом объединения или пересечения фигур (рис. 1):

>> mx=1; my=mx*2/3; sigX=3; sigY=sigX*0.8; r=0.6; a=5; b=a*0.7;

>> X=Norm_2([mx;my],[sigX,sigY],r);

>> R=Rect([a,b]); C=Circ(b/2, [0;b/2]);

>> pC=Ver(X,C), pRC = Ver(X, R&C)

>> ShowAll(C,R,R&C,'Fc',X)

pC =

    0.1684

pRC =

    0.1069

Рис. 1. Иллюстрация операций с объектами фигур и вычисления вероятностей попадания

Результат выполнения электронной формулы Ver зависит от структуры объекта G. Если G – массив фигур, P – вероятность попадания хотя бы в одну из них, p – массив вероятностей попадания в отдельные фигуры. Из массива вероятностей p конструктор класса случайных событий Randevent создает массив событий с вероятностями из массива p: A = Randevent(p). Выражение для суммы событий sum(A) в данном примере даст результат, совпадающий со значением P, так как элементы массива A – несовместные события.

Технология электронных формул для вероятностных расчетов

Вычисление вероятности сложных событий выполняют объекты класса событий Randevent переопределенными аддитивными и мультипликативными операциями, например:

>> [A,B,C]= Randevent(0.3, 0.4, 0.5); D=A+B*C, E=(A+B)*C

Событие D: P(D) =   0.44

Событие E: P(E) =   0.29

В операциях умножения учитывается зависимость между событиями, заданная условными вероятностями P(A|C) = 0.45, P(B|C) = 0.35:

>> A=Set(A,C,0.45); B=Set(B,C,0.35); D=A+B*C, E=(A+B)*C

Событие D: P(D) =  0.396

Событие E: P(E) =  0.321

В операциях с массивами событий функция sum (сложение) учитывает совместность слагаемых, функция prod (произведение) – попарную зависимость:

>> p=[0.2 0.4 0.45 0.25 0.5 0.1 0.6 0.15 0.4];

>> SYS = Randevent(p), A = sum(SYS), B=prod(SYS(1:4))+prod(SYS(5:9))

Группа событий SYS: P(SYS(1:9)) =    0.2   0.4  0.45  0.25   0.5   0.1   0.6  0.15   0.4

Событие A: P(A) =  0.982

Событие B: P(B) = 0.0108

Электронная формула Bayes класса Randevent вычисляет апостериорную вероятность интересующей гипотезы по последовательности событий (A, A, …., not A, …), наблюдавшихся в контрольных испытаниях. В качестве первого аргумента она получает объект уточняемой гипотезы Hi и последовательность объектов A и not(A), возвращаемый результат – апостериорная вероятность уточняемой гипотезы:

HiPost=Bayes(Hi, A, …, not(A), …).

Порядок чередования событий и их отрицаний в списке аргументов безразличен, можно заменить группу одинаковых событий парой аргументов – (событие, кратность):

HiPost=Bayes(Hi, A, n, not(A), m).

Вычисления по частной, общей и обобщенной формулам Бернулли выполняет электронная формула Ber:

out = Ber(p, n, m),

где p – вероятность успеха в каждом испытании, или вектор-строка вероятностей успеха или вектор-столбец вероятностей каждого из множества возможных исходов; n – число испытаний; m – интересующее число успехов, или вектор-строка диапазона интересующего числа успехов или вектор-столбец интересующей комбинации чисел исходов. Выходная величина out в прямых задачах имеет значение вероятности наступления m событий, в обратных задачах – принимает значение вычисленной величины, заданной пустым вектором [ ] в списке аргументов. Электронная формула Ber с пустым первым аргументом и заданной четвертым аргументом доверительной вероятностью Dov вычисляет доверительный интервал [p1, p2] для вероятности успеха по m успешным исходам в n испытаниях.

Дискретные случайные величины

В объектно-ориентированной вероятностной модели объекты распределения участвуют в алгебраических операциях, сами вычисляют МО и СКО, вероятности событий. Общие свойства дискретных распределений (ряд распределения и числовые характеристики) инкапсулированы в базовый класс DRV дискретных СВ. Теоретические законы распределения, различающиеся набором параметров и их интерпретацией, реализованы в производных классах от DRV.

Вызов конструктора DRV без аргументов заставляет его вывести сведения о себе и основных функциях базового класса, а также список производных классов дискретных СВ:

>> DRV

 конструкторы, создающие объект дискретного распределения X:

  X=DRV(x;p) или X=DRV([x;p])-  x - массив значений; p - массив вероятностей

  X=DRV(A); A-случайное событие (объект класса Randevent), X-индикатор события A

 основные функции класса:

  Val(X,Ind) - возможные значения СВ X (с индексами Ind, если 2-й аргумент задан)

  Ver(X,Ind) - вероятности значений СВ X (с индексами Ind, если 2-й аргумент задан)

  MO(X), SKO(X) - вывод МО и СКО СВ X)

  Show(X1,s1,X2,s2,...) - вывод многоугольников распределения СВ X1,X2,...)

 конструкторы производных классов:

  BIN(p,n) - биномиальное распределение

  POI( L ) - распределение Пуассона

  GEO(p,d) - геометрическое распределение (d=0) или сдвинутое геометрическое (d = 1, 2, …)

  HyperGEO(N,R,L) - гипергеометрическое (R дефектных из N, L - объем выборки)

Массив вероятностей, принимаемых конструктором, проверяется на нормированность единицей. При этом, если один из элементов массива вероятностей задан числом -1, он заменяется дополнением суммы остальных элементов до единицы. Если модуль отрицательного числа меньше единицы, элемент заменяется абсолютным значением, а дополнение до единицы равномерно распределяется между нулевыми элементами массива. Аргументом конструктора DRV может быть случайное событие, тогда объект создается как индикатор этого события.

Алгебраические операции с объектами СВ выполняются по правилам теории вероятностей, сохраняют свойства распределений, такие как устойчивость закона Пуассона к сложению:

>> X2=POI(2); X3=POI(3); X5=X2+X3

X5: Распределение Пуассона (5) MO=5, SKO=2.24

Базовый класс непрерывных случайных величин

Базовый класс CRV непрерывных СВ содержит закон распределения в виде функции плотности вероятности или функции распределения, или их дискретных значений на сетке. Объекты непрерывных СВ создает конструктор CRV:

>> CRV

 конструкторы непрерывного распределения, создающие объект X:

  X=CRV(x;p) -  x - массив значений; p - массив вероятностей)

  X=CRV(S;lim) -  S – формула F(x) или f(x); lim – [1x2]-вектор границ интервала значений

 конструкторы классов для основных законов распределения:

  EXP( L )    - показательное распределение

  ERL( L,k ) - распределение Эрланга порядка k

  RND(a,b) - равномерное распределение

  RAYL(p)   - распределение Рэлея

  WEI(p1,p2) - распределение Вейбулла

  NORM(m,s)- нормальное распределение

  GAMMA(a,b) - Гамма-распределение

Базовый класс CRV создает объекты произвольных распределений, заданных формулой или статистическим рядом, объекты распределения строят графики интегрального и дифференциального законов:

>> XF=CRV('2/pi*asin(x/10)',[0,10]), Xf=CRV('2/pi./sqrt(100-x.^2)',[0,10]),Show(XF,'F:b',XF,'f:r')

XF: Распределение F:2/pi*asin(x/10) [0 10]: MO=6.36, SKO=3.07

Xf: Распределение f:2/pi./sqrt(100-x.^2) [0 10]: MO=6.36, SKO=3.08

В качестве своего значения объекты выводят тип или происхождение, интервал возможных значений, математическое ожидание и среднеквадратичное отклонение.

Показательный закон и связанные с ним классы распределений

Вероятностный анализ событий, связанных с пуассоновскими потоками отказов, заявок, сигналов может быть выполнен без методических упрощений объектами класса показательного закона EXP и родственных классов закона Эрланга ERL и Вейбулла WEI. В следующем примере конструктор EXP создает объект показательного распределения с параметром 2, объекты выполняют алгебраические действия по правилам теории вероятностей, функция Show строит графики плотности и функции распределения:

>> T=EXP(1.5), Show(T,'Fk', T, 'fr'), T3=T+T+T, figure, Show(T,'k',T+T,'r',T3,'b',T3+T,'g')

T: Распределение экспоненциальное (1.5)  [0 Inf]: MO=0.67, SKO=0.67

T3: Распределение Эрланга порядка 2(1.5) MO=2, SKO=1.15

Происхождение классов непрерывных распределений от базового класса сеточного распределения обеспечивает результативность операций с объектами. Массив объектов X класса EXP, созданный массивом параметров, аппроксимирует произведение Prodf(X) = Y объектом класса WEI, выводит графики функций распределения объектов X, плотности и функции распределения СВ Y, а также функцию распределения Вейбулла (точками), близкую к графику функции распределения СВ Y:

>> X=EXP(1.2:0.2:1.8), Y=Prodf(X), W=WEI(Y), Show(X,'Fb', Y, 'r',Y,'Fk',W,'Fr.')

Массив X:

 Распределение экспоненциальное (1.2)  [0 Inf]: MO=0.83, SKO=0.83

 Распределение экспоненциальное (1.4)  [0 Inf]: MO=0.71, SKO=0.71

 Распределение экспоненциальное (1.6)  [0 Inf]: MO=0.63, SKO=0.63

 Распределение экспоненциальное (1.8)  [0 Inf]: MO=0.56, SKO=0.56

Y: Распределение  [0, Inf]: МО =1.37, SKO = 0.83

W: Распределение Вейбулла (0.6416,1.7807) MO=1.39, SKO=0.81

Точность аппроксимации можно оценить сравнением с вероятностью отказа или безотказной работы в заданном интервале, вычисленной по формуле сложного события в классе Randevent.

Нормальный закон распределения

Сведения об основных функциях класса NORM можно получить по имени класса.

>> NORM

 X=NORM(m,s) - конструктор объектов нормального распределения

 X=NORM(m) или X=NORM - конструктор объектов (по умолчанию m=0, s=1)

   основные функции класса:

 y = f(X,x) - вычисляет плотность y распределения СВ X на сетке x

 x = Net(X,n,N) - строит сетку n точек на интервале "N сигм" (по умолчанию n=50,N=4)

 T = Gen(X,m,n,Y) - генерирует nxm-матрицу T реализаций X, смещенных на Y (n=1,Y=[])

 F=Ver(X, X<a) - вероятности события (X<a), Ver(X,X>a) – вероятность события (X>a)

 F=Ver(X, [a,b]) - вероятность попадания в интервал

 F=Ver(X, str) - вычисление по полной вероятности str - выражение условной вероятности

Класс двумерных нормально распределенных векторов

Конструктор Norm_2 создает объект, выполняющий все действия, предусмотренные математической моделью двумерных нормальных распределений. Класс Norm_2 наследует базовый класс геометрических объектов вместе с методами пространственных преобразований (сдвиг, масштабирование, повороты). Вектор МО как центр рассеивания позиционирует объект, СКО определяет размеры единичного эллипса рассеивания. Двумерный вектор-столбец принимается конструктором как вектор МО, положительный [1, 2]-вектор – как СКО. Скалярный аргумент учитывается конструктором как коэффициент корреляции. Сведения о способах определения объектов и методах работы с ними можно получить вызовом конструктора Norm_2 без аргументов и принимающей переменной:

>> Norm_2

   Конструкторы объектов нормального распределения:

X=Norm_2(m, K); m - 2x1-вектор МО =[0;0], K – корреляционная матрица

X=Norm_2(m, s, r);  m=[0;0], s - 1x2-вектор СКО =[1;1], r - коэф. корреляции =0

X=Norm_2(X1, X2, r); X1, X2 - объекты класса NORM, r = 0

   Основные функции класса:

 dens = f(X, x, y) - вычисляет плотность dens распределения СВ X на сетке x, y

 [x, y] = Net(X,n,N) - строит сетку n точек на интервале "N сигм" (по умолчанию n=50,N=4)

T = Gen(X,m,n) - генерирует nxm-матрицу T реализаций X    

Y = Rot(X,a) - поворот осей на угол a1[град]

P = Ver( X, G) – вычисляет вероятность P попадания случайной точки X в область G 

Переменную в левой части конструктор без параметров определяет объектом стандартного распределения с единичными СКО без систематических ошибок и корреляции. Операции сдвига, масштабирования и поворота иллюстрируют операции с объектом Y0:

>> Y0=Norm_2, Y1 = Y0*[2, 1], Y2 = Y1+[12;3], Y3 = Rot(Y1,-30) + [12; 3]

 Norm_2 Y0: MO = [0 0], CKO = [1 1]

 Norm_2 Y1: MO = [0 0], CKO = [2 1]

 Norm_2 Y2: MO = [12 3], CKO = [2 1]

 Norm_2 Y3: MO = [12 3], CKO = [1.8028 1.3229], r = 0.54

Электронная формула Ver(X, G) вычисляет вероятность попадания нормально распределенной точки X в геометрическую область G численным интегрированием или по специальным формулам в частных случаях.

Функции случайных величин

Аналитические преобразования функции случайных величин могут быть заменены универсальной процедурой, обрабатывающей интерпретируемое выражение функции согласно теоретической формуле. Электронная формула FUN получает выражение функции S, один или несколько случайных аргументов X1, ..., параметры a1, ..., и создает объект распределения Y класса CRV: Y=FUN(S, X1, ..., a1, ...)

Выражение S для площади проекции параллелепипеда с параметрами a, b, c содержит перечень параметров в том порядке, в каком они принимаются из списка входных данных, угол Fi, равномерно распределенный определен объектом класса равномерного распределения RND, объект распределения угла Teta – по закону косинуса:

>> S = 'b*c*sin(Teta) + a*(c*abs(sin(Fi)) + b*abs(cos(Fi))).*cos(Teta) : a,b,c'; a=9; b=6; c=4;

>> Teta = CRV('cos(x)', [0, pi/2]); Fi = RND([0, 2*pi]); Pr3 = FUN(S, Teta, Fi, a, b, c)

Pr3: Распределение  [24.0105, 69.1915]: МО =56.96, SKO = 9.42

Для вычисления площади проекции S геометрического объекта G класса многогранников PolyGran на плоскость, ориентированную под сферическими углами Teta, Fi, можно использовать электронную формулу ProectS: S = ProectS(G, Teta, Fi)

Мультипликативные функции случайных величин

Мультипликативные операции с объектами распределений учитывают особенности сомножителей. В примере построенные законы распределения площади прямоугольника, длины сторон которого X1 ∈ N(10, 2) и X2 ∈ N(15, 4) в одном случае независимы, а двух других зависимы с коэффициентами корреляции 0,5 и  –0.5:

>> r =    0; X=Norm_2([10;15],[2 4], r); [X1,X2]=X12(X); Y1=X1*X2; Show(Y1, 'b')

>> r =  0.5; X=Norm_2([10;15],[2 4], r); [X1,X2]=X12(X); Y2=X1*X2; Show(Y2, 'r')

>> r = -0.5; X=Norm_2([10;15],[2 4], r); [X1,X2]=X12(X); Y3=X1*X2; Show(Y3, 'k')

Аддитивные функции случайных величин

Сумма двух независимых случайных величин (композиция), равномерно распределенных в одном интервале, подчиняется закону Симпсона, сумма нескольких таких величин имеет нелинейную плотность распределения, приближающуюся к нормальному закону при числе слагаемых более пяти:

>> R1=RND(0,1),R2=R1+R1, R3=R2+R1,R4=R3+R1, R6=R4+R2, N=NORM(R6), Show(R1,R2,R3,R4,R6,N ,'.')

R1: Распределение RND [0 1]: MO=0.5, SKO=0.29

R2: Распределение Симпсона[0, 2]: МО =1, SKO = 0.41

R3: Распределение  [0, 3]: МО =1.5, SKO = 0.5

R4: Распределение  [0, 4]: МО =2, SKO = 0.58

R6: Распределение  [0, 6]: МО =3, SKO = 0.71

N: Распределение нормальное (3 0.71) : MO=3, SKO=0.71

Композиция нормального и равномерного законов выполняется операцией сложения с объектами классов NORM и RND:

P = Ver(NORM(m, sigma) + RND(a,b), [z1, z2])

Применение классов в теории массового обслуживания

Универсальную модель занятости многоканальной СМО с отказами и очередями воспроизводит стохастическая процедура, реализованная в программе Qsystem. Она получает объекты потоков заявок X, обслуживания Y и число каналов n. Признаком СМО с ожиданием служит отрицательный знак аргумента n. Дополнительными необязательными аргументами задаются число повторений процесса N для набора статистики (по умолчанию, 10000), длительность процесса T (10) и число net (50) моментов времени, в которые оценивается вероятность наличия свободного канала. Программа N раз имитирует последовательность поступления заявок, длительность их обслуживания, отказы в обслуживании или постановку в очередь. По накопленной статистике восстанавливается зависимость незанятого состояния СМО с отказами или среднего времени ожидания в СМО с очередями. Для аппроксимации зависимости подбирается подходящая функция, которая выводится в качестве результата вместе с параметрами.

Вычисление вероятностей попадания в эллипсоид по универсальной электронной формуле проще, чем по специальной формуле, использующей таблицы функции Лапласа:

>> k =[0.5;1;2;3;4]; X=CHI2(3); p=Ver(X,k.^2); g = 2*FLaplas(k)-sqrt(2/pi).*k.*exp(-k.^2/2); p=p',g=g'

p =   0.030857      0.19873      0.73851       0.9707      0.99887

g =   0.030857      0.19873      0.73851       0.9707      0.99887

Электронная формула Ver(X, G) заменяет интегрирование по области G суммированием элементов вероятностей в малых окрестностях N точек r_i, равномерно распределенных в G. Расчетные точки в заданном количестве N равномерно распределяет в области G сам объект геометрии, он же вычисляет площадь S_D = Area(G), объект распределения X вычисляет в этих точках значения плотностей распределения:

Pnt = Gen(X, N); P = sum( Pdf(X, Pnt) ) * Area(G) / N

Расчет надежности сложных систем

Расчет надежности сложных схем без упрощений и оптимизацию надежности легко выполнить методами классов непрерывных случайных величин (СВ). Так как все "именные" классы распределений (EXP, ERL, WEI, и др.) имеют общий базовый класс CRV с сеточной функцией распределения, математические действия с объектами СВ результативны. В частности, перемножение функций распределения выполняет функция Prodf, создающая объект класса CRV. Следующая команда определяет массив объектов X класса показательного распределения EXP массивом параметров, аппроксимирует произведение Y = Prodf(X) объектом класса WEI (Вейбулла), выводит графики функций распределения объектов X, плотности и функции распределения СВ Y, а также функцию распределения Вейбулла (точками), близкую к графику функции распределения СВ Y (рис. 2):

>> X=EXP(1.2:0.2:1.8), Y=Prodf(X), W=WEI(Y)

>> Show(X,'Fb', Y,'Fk',W,'Fr.')

Массив X:

 Распределение экспоненциальное (1.2)  [0 Inf]: MO=0.83, SKO=0.83

 Распределение экспоненциальное (1.4)  [0 Inf]: MO=0.71, SKO=0.71

 Распределение экспоненциальное (1.6)  [0 Inf]: MO=0.63, SKO=0.63

 Распределение экспоненциальное (1.8)  [0 Inf]: MO=0.56, SKO=0.56

Y: Распределение  [0, Inf]: МО =1.37, SKO = 0.83

W: Распределение Вейбулла (0.6416,1.7807) MO=1.39, SKO=0.81

Рис. 2. Функции отказов элементов,

параллельного звена и аппроксимация законом Вейбулла (точки)

Вероятность безотказной работы параллельной схемы в течение 2-х ед. времени можно вычислить по точной формуле P(T > 2) = 1 – ∏F_i(2) и с помощью виртуальной функции Ver, применяемой к объектам Y и W:

>> p2 =1- prod(1-exp(-[1.2:0.2:1.6]*2)), pY = Ver(Y,Y>2), pW=Ver(W,W>2)

p2 = 0.1808  pY = 0.2033 pW = 0.2103

Аппроксимация функции отказов параллельного звена законом Вейбулла (объект W) не дает вычислительных преимуществ по сравнению с объектом Y и уступает в точности, но позволяет заменить последовательно-параллельную схему соединения последовательной. Так, последовательно-параллельную схему (рис. 3, а),  логика отказов которой определена текстовым выражением S, а параметры – массивом Lam, можно заменить эквивалентной последовательной схемой и построить график функции распределения отказов:

>> S='1 or 2 or 3 and 4 or 5 and 6 and 7 or 8 or 9 and 10';  Lam = [3, 2, 12, 15, 10, 11, 11, 1, 18, 20]/100;

>> E=EXP(Lam); W34=WEI(Prodf(E(3:4))); W57=WEI(Prodf (E(5:7))); W910=WEI(Prodf (E(9:10)),'F');

С помощью функции MyPar(X, k), которая возвращает параметр масштаба (k = 1) или параметр формы (k = 2), выделим параметры закона Вейбулла λ и α из эквивалентных законов для резервированных звеньев W34, W57, W910, сформируем векторы параметров эквивалентных элементов:

>> L=[Lam(1),Lam(2), MyPar(W34,1) , MyPar(W57,1) ,Lam(8), MyPar(W910,1)];

>> A=[1,1, MyPar(W34,2) , MyPar(W57,2) ,1, MyPar(W910,2)];

Функцию распределения эквивалентной последовательной схемы вычислим на расчетной сетке в интервале [0, 20] (рис. 3, б – сплошная линия):

>> T=0:0.1:20; Z=[]; for t=T Z(end+1)=1-exp(-sum((L*t).^A)); end; plot(T,Z), F=CRV(Z,T)

F: Распределение CRV [0 0.999]: MO=0.71, SKO=0.24

Рис. 3. Последовательно-параллельное соединение элементов системы (а) и функция распределения отказов эквивалентной последовательной схемы (б)

Сравним вероятность наступления отказа в течение 1-й ед. времени по упрощенной схеме с точной вероятностью, вычисленной по формуле сложного события в классе Randevent согласно логической схеме:

>> P1=1-exp(-sum((L*1).^A)), B=Randevent(DF(E, 1)); R1=Set(B, S)

P1 =  0.1454

Событие R1: P(R1) =  0.102

Существенное отличие вероятности отказа 0.1454 по упрощенной схеме от истинной 0.102 объясняется погрешностями аппроксимации надежности параллельных звеньев законом Вейбулла. Точную функцию отказов можно построить, вычислив истинные вероятности отказа в интервалах разной длительности в классе Randevent (рис. 3, б – пунктирная линия):

>> Ft=[];T=0:0.1:20;for t=T; B=Randevent(DF(E, t)); Ft(end+1)=Value(Set(B,S)); end, hold on, plot(T,Ft,'r--')

Функция распределения Ft на сетке T, свободная от погрешностей аппроксимации, адекватно представляет резервированное звено одним элементом, эквивалентным по надежности:

>> Elem = WEI(T,Ft), Show(Elem,'F:k.')

Elem: Распределение Вейбулла (0.2147,1.4501) MO=4.22, SKO=2.96

Корректное упрощение последовательно-параллельной схемы позволяет распространить объектное моделирование надежности на системы более сложной структуры с сохранением высокой точности прямого вычисления надежности.

Достижение требуемого уровня надежности проектируемой системы связано с обоснованием необходимой надежности элементов, рациональным выбором структуры системы, оптимальным резервированием. Решение этих проблем, как и все обратные задачи, требует адаптации расчетных моделей, способствующей применению эффективных поисковых методов.

Методы оптимизации надежности реализованы в классе надежности Reliab, который не входит в иерархию классов распределений, но использует классы EXP, WEI для выполнения прямых вычислений. Структура системы и параметры надежности элементов передаются конструктору Reliab одним массивом ячеек так, что каждая группа последовательных (параллельных) соединений помещается в отдельный вложенный массив ячеек в виде строки (столбца):

>> R=Reliab(1,{0.03 0.02, {0.12; 0.15}, {0.1; 0.11; 0.11}, 0.01, {0.18; 0.2}}), MyPar(R)

Надежность R = 0.8983 (1)

  0.0300    0.0200    0.1200    0.1500    0.1000    0.1100    0.1100    0.0100    0.1800    0.2000

  1 or 2 or 3 and 4 or 5 and 6 and 7 or 8 or 9 and 10

Массив ячеек для создания объекта R описывает приведенную выше схему, поэтому две последних строки, сформированные командой MyPar(R), совпали с исходными массивами S и Lam для этой схемы, а надежность R = 0,8983 – с дополнением к вероятности отказа W=0,102. Конструктор Reliab получил также продолжительность рабочего периода (1) для вычисления вероятностей отказа. В качестве своего значения объекты класса выводят дополнение к вероятности отказа в рабочий период – надежность. Элементами исходного массива ячеек могут быть как параметры потока отказов, так и вероятности безотказной работы. При определении объекта массивом вероятностей продолжительность рабочего периода не используется, скалярный аргумент должен отсутствовать или иметь неопределенное значение (Inf):

>> P=Reliab( 0,{0.68 0.72, 0.87 0.62 0.78, 0.65})

Надежность P = 0.1339 (Inf)

Логические схемы любой сложности одинаково просто представлять структурированными выражениями в массиве ячеек. Роль постановщика задачи сводится к формированию объектов и структурированию исходных данных, расчеты выполняются взаимодействием объектов, без предварительной алгоритмизации.

Пример семинара

Случайные события

1.1. Непосредственный расчет вероятностей

Задача 1. На отрезок единичной длины случайно брошены две точки. Какова вероятность того, что расстояние между точками меньше 0,25?

Проверка:

>> Om = Rect([0;0],[1;1]),A=Polygon2([0,0,0.75,1,1,0.25,0;0,0.25,1,1,0.75,0,0]);

>> P=Area(A)/Area(Om), Show(Om,A,'Fr'), daspect([1,1,1])

P = 0.4375

Задача 2. Плоская фигура составлена из прямоугольников 6×1.5 с центром в точке (0,1.25) и 2×1 с центром в точке (0,2.5), и двух кругов диметром 1 с центрами в точках (-1.5,0) и (1.5, 0). Случайная точка равномерно распределена в границах фигуры. Вычислить вероятность попадания точки в меньший прямоугольник или круг.

Проверка:

>> K=Rect([-3;0.5],[3;2]),B=Rect([-1;2],[1;3]),C1=Circ(0.5,[-2;0.5]);C2=C1+[4;0];

>> Om=K|B|C1|C2; P=(2*Area(C1)+Area(B))/Area(Om)

>> Show(K,'r',B,C1,C2), daspect([1,1,1])

P = 0.3030

1.2. Вероятности сложных событий

Задача 3. События A, B, С независимы, Р(А) = 0.3, Р(B) = 0.4, Р(C) = 0.5. Вычислить вероятность событий D = A + ВC и E = (A + В)C.

Проверка:

>> [A,B,C]=Randevent(0.3,0.4,0.5); D=A+B*C, E=(A+B)*C

Событие D: P(D) = 0.44

Событие E: P(E) = 0.29

Задача 4. Вероятности событий A, B, С те же, что и в предыдущей задаче, но теперь события B, С зависимы. Известна условная вероятность Р(В|C) = 0.45. Как изменятся вероятности событий D и E?

Проверка:

>> [A,B,C]=Randevent(0.3,0.4,0.5); B=Set(B,C,0.45); D=A+B*C, E=(A+B)*C

Событие D: P(D) = 0.458

Событие E: P(E) = 0.308

1.3. Формула Байеса в принятии решений

Задача 5. Событие A – обнаружение дефекта в случайно выбранном изделии, гипотеза H₁ определяет характерную для данного производства надежность (отсутствие дефектов) с вероятностью P(H₁) = 0,8. Если событие A наступило (обнаружен дефект в случайно выбранном изделии), условная вероятность гипотезы P(H₁|A) может стать меньше критической 0,95. Методы контроля качества продукции могут ошибочно забраковать кондиционное изделие с вероятностью P(A/H₁) = 0,1, а вероятность обнаружить имеющийся дефект (гипотеза H₀) составляет P(A/H₀) = 0,8. Спланировать процедуру выборочного контроля для определения кондиционности партии.

Указания:

Определив гипотезы и условные вероятности события в классе Randevent, получим результат с помощью электронной формулы Bayes:

>> [H0,H1]=Randevent(0.2,0.8), A=SetDepend(H0,0.8,H1,0.1); H=Bayes(H1, not(A))

Событие H: P(H) =  0.947

Если дефект был обнаружен, годность партии маловероятна:

>> H=Bayes(H1, A)

Событие H: P(H) =  0.333

Не исключено, однако, подтверждение годности после дополнительного количества успешных испытаний. В формулу Bayes нужно подставлять апостериорную версию гипотезы с событиями дополнительных испытаний, либо априорную гипотезу с полной последовательностью событий (что предпочтительнее, так как объект события A связан с H1):

>> H=Bayes(H1, A, not(A), not(A))

Событие H: P(h) =   0.91

Результат двух успешных испытаний после одного неудачного оказался ниже доверительной вероятности, для подтверждения годности партии необходимо продолжить испытания. Использование исходной гипотезы H1 с априорной вероятностью предпочтительнее, так как последовательность аргументов соответствует последовательности наблюдаемых событий:

>> H=Bayes(H1, A, not(A), not(A), not(A))

Событие H: P(H) =  0.979

Последовательность испытаний в данном случае подтверждает кондиционность партии. Порядок чередования событий и их отрицаний в списке аргументов безразличен, можно группировать одинаковые аргументы и заменить группу парой аргументов – событие, кратность:

>> H_13=Bayes(H1, A, not(A),3), H_10=Bayes(H1, A), H_24=Bayes(H1, A, 2, not(A),4)

Событие H_13: P(H_13) =  0.979

Событие H_10: P(H_10) =  0.333

Событие H_24: P(H_24) =  0.962

Согласно методике контроля, испытания могут проводиться в два или более этапов с определенным числом повторений. Если на первом этапе гипотеза H1 не принимается, например:

>> H=Bayes(H1, A, 1, not(A), 2)

Событие H: P(H) =   0.91

проводится дополнительная серия, по результатам которой партия принимается или бракуется:

>> H=Bayes(H1, A, 1+1, not(A), 2+2)

Событие H: P(H) =  0.962

1.4. Случайные величины

Задача 6. Вычислить средний расход снарядов в стрельбе до первого попадания при вероятности попадания в одном выстреле p₁=0,2. Число снарядов в боекомплекте n=12; n=50.

Решение:

Построим сдвинутое на 1 геометрическое распределение и применим функцию min для возможного расхода:

>> p=0.2;n=12; X=GEO(p,1); X12=min(X,12),X50=min(X,50)

X12: Распределение дискретное MO=4.66, SKO=3.46

X50: Распределение дискретное MO=5, SKO=4.47

Задача 7. Разыграть 50000 случайных реализаций длины проекции стержня в интервале [0, l] и построить две гистограммы с шириной разрядов 0.5 и 0.25. На том же графике показать функцию плотности распределения проекции.

Решение:

>> L=10; N=50000; x=rand(1,N)*pi/2;X=sin(x)*L; [F,f,H]=SmartHist(X,[1,10],20); Show(H,'r',f,'r'), hold on, [F,f,H]=SmartHist(X,[1,10],40); Show(H,f)