Software & Systems

Известно, что марковские (полумарковские) модели [1] являются удобным инструментальным средством для исследования сложных, в том числе производственно-экономических, систем [2, 3]. Указанные модели подобных систем обычно представляются в виде ориентированного графа, вершины которого соответствуют возможным состояниям системы, а веса соединяющих их дуг – таким числовым параметрам, как интенсивности переходов из одного состояния в другое (интенсивностям потоков событий). При моделировании в качестве показателя эффективности обычно применяется значение вероятности нахождения системы в некотором состоянии, при этом основным затруднением в использовании данных моделей является неполнота статистической информации о значениях интенсивностей [3]. В статье рассмотрена задача нахождения оценки диапазона возможных значений вероятности пребывания системы в заданном состоянии при интервальной неопределенности об интенсивностях потоков событий.

Постановка задачи. Предполагается, что в рамках марковской модели задан граф состояний системы со стационарными пуассоновскими потоками событий, переводящими систему из состояния si в состояние sj. Число вершин графа ограничено, n£10. Точные значения интенсивностей переходов wij неизвестны (очевидно, все wij³0), для ненулевых wij заданы лишь интервалы их возможных значений wij min£wij£wij max, при этом предполагается, что длина каждого интервала намного меньше, чем положение его центра, то есть

wij max–wij min<<(wij max+wij min)/2,

wij max>0; i=1,2,…, n; j=1,2,…, n.                       (1)

Заметим, что общее количество N ненулевых параметров wij в случае, если каждая из n вершин графа соединена дугами с остальными (n-1) вершинами, равно N=n×(n-1); на практике обычно имеют место 1–3 связи каждой вершины с другими, так что можно полагать N»2×n.

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

Метод решения задачи. При известных значениях интенсивностей wij, то есть в условиях полной информации о параметрах модели, вероятности переходов из одного состояния системы в другое в установившемся режиме (предельные вероятности) находятся из уравнений Колмого- рова [4]:

 (i = 1,2,…, n).                 (2)

Первая сумма в левой части формулы (2) распространяется на те значения j, для которых возможен непосредственный переход из состояния sj в состояние si (то есть для которых wji¹0), а вторая – на те значения j, для которых возможен непосредственный переход из si в sj (то есть wij¹0). 

Приведенные уравнения дополняются нормировочным условием

.                                                                   (3)

Совокупность уравнений (2) и (3) позволяет определить значения вероятностей p1, p2, …, pn, при этом для получения однозначного решения нормировочное условие (3) можно использовать вместо любого из уравнений (2), например, вместо последнего n-го уравнения. Совокупность используемых уравнений можно при этом записать в векторно-матричной форме:

,                                                                  (4)

где  – вектор-столбец, элементами которого являются искомые вероятности (здесь T – символ транспонирования); Аn´n – матрица, элементы первых (n-1) строк которой определяются в соответствии с выражениями (2) и зависят, следовательно, от значений интенсивностей wij, а элементы n-й строки – единицы;  – вектор-столбец, все элементы которого, кроме последнего, нулевые, а последний, n-й – единица.

Отсюда получаем:

.                                                               (5)

Если нас интересует вероятность pi только какого-то одного i-го состояния, то на основании (5) можно записать:

,                                                         (6)

где  – вектор-строка, у которого все элементы, кроме i-го, нулевые, а ci=1.

На практике решение указанной задачи затруднено вследствие недостатка статистических данных, на основе которых можно было бы определить вероятностные характеристики марковской модели. В условиях рассматриваемой задачи имеющаяся неопределенность о параметрах модели отражена совокупностью неравенств (1), задающих интервалы возможных значений этих параметров.

Точное решение поставленной задачи достигается решением двух оптимизационных задач (нелинейного программирования):

,                               (7)

,                               (8)

где {wij} – совокупность параметров модели (интенсивностей); Ww – область допустимых значений этих параметров (гиперпараллелепипед, заданный совокупностью неравенств wij min£wij£ £wij max), а обозначение A({wij}) указывает на зависимость матрицы A от искомых параметров.

Нетрудно показать, что аналитическое решение (7), (8) удается получить лишь в простейших случаях для небольших размерностей модели, n£3. Однако и численное решение, например, с помощью современных систем компьютерной математики типа Mathcad, MATLAB, возможно лишь при числе аргументов (ненулевых параметров wij) не более 12–15.

При большей размерности задачи, по-видимому, единственное, что можно сделать, – это получить лишь некоторую оценку интервала возможных значений pi.

Для получения такой оценки введем обозначения:

– для средней точки i,j-го ненулевого параметра модели

;                                            (9)

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

.                       (10)

Отметим, что, исходя из соотношений (2) и (3) для элементов матрицы A, имеем:

A({})=A({})+A({})=+DA,   (11)

где =A({}), DA=A({}).

Обозначим через  решение системы

A({})×=×=,                                         (12)

а через x+Dx – решение системы

A({})×(+Dx)=

=(+DA)×(+Dx)=.                                        (13)

Из последнего равенства будем иметь

×+DA×+×Dx+DA×Dx=.                      (14)

Вычитая из этого соотношения (12), получим

×Dx = – DA× – DA×Dx,                                   (15)

отсюда 

Dx=DA× + DA ×Dx).                              (16)

Теперь, используя норму 1 матриц и векторов и на основании свойств нормы [5], запишем:

).

Норма 1 вектора  есть сумма абсолютных величин его элементов. Поскольку эти элементы – неотрицательные величины, являющиеся вероятностями полной группы событий, то  (см. соотношение (3)). Ввиду этого далее следует

(1 –   

и, если дополнительно предположить, что  < 1, .

Умножим обе части полученного неравенства на норму вектора-строки  (см. выше):

.

Поскольку

и на основании определения вектора , , окончательно имеем

                                  (17)

(знак модуля в правой части неравенства опущен, поскольку левая часть всегда неотрицательная).

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

pi min=pi–Dpi,                                                            (18)

pi max=pi+Dpi,                                                     (19)

где , а  находится в результате решения уравнения (12).

Иллюстрирующий пример. Пусть задана структура вида, представленного на рисунке, необходимо найти предельную вероятность p2 пребывания системы в состоянии s2. 

Пусть заданы следующие интервалы для параметров системы:

1.975 £ w12 £ 2.025,

0.975 £ w14 £ 1.025,

0.975 £ w24 £ 1.025,

2.950 £ w31 £ 3.050,

1.975 £ w41 £ 2.025,

1.975 £ w43 £ 2.025.

Тогда в соответствии с изложенным имеем:

A,          (20)

=[0 0 0 1]T, =[0 1 0 0], =2, =1, =1, =3, =2, =2;

=,

DA=.

Использование уравнения (12) и соотношения  дает p2=0.471, а из неравенства (17) следует Dp2£0.118 (расчеты проводились в среде MATLAB версии R2008a) и далее – по (18), (19) – получаем окончательно p2 min³0.353, p2 max£0.589.

Проверка полученных результатов проводилась путем нахождения решений оптимизационных задач (7), (8) с помощью функции MATLAB fmincon (функции минимизации скалярной функции многих аргументов [6]). В этом случае интервал возможных значений вероятности p2 таков: 0.459£p2£0.482.

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

В заключение отметим, что рассмотрена марковская модель производственно-экономической системы, отличающаяся интервальной неопределенностью об интенсивностях потоков событий. Целью исследования являлось нахождение диапазона возможных значений вероятности пребывания системы в заданном состоянии. Показано, что такая задача сводится к задаче нелинейной оптимизации скалярной функции многих переменных. Указано на трудность получения как аналитического решения данной задачи, так и прямого использования численных методов ввиду возможной большой размерности задачи. Предложен метод приближенного нахождения границ интервала неопределенности для требуемой вероятности, использующий аппарат матричной алгебры, который с вычислительной точки зрения представляется более простым, чем прямое решение оптимизационной задачи. Действительно, операции даже с матрицами размера 10´10 проще, чем поиск экстремума функции 20 и более переменных при наличии ограничений типа неравенств.

Литература

1. Гнеденко Б.В., Коваленко И.Н. Введение в теорию массового обслуживания. М.: Наука, 1987.

2. Куликов Г.Г., Флеминг П.Дж., Брейкин Т.В., Арьков В.Ю. Марковские модели сложных динамических систем: идентификация, моделирование и контроль состояния (на примере цифровой САУ ГТД). Уфа: УГАТУ, 1998.

3. Бояринов Ю.Г., Мищенко В.И. Основные направления повышения эффективности полумарковских моделей производственно-экономических систем // Программные продукты и системы. 2009. № 2. С. 144–148.

4. Вентцель Е.С., Овчаров Л.А. Теория случайных процессов и ее инженерные приложения. М.: Наука, 1991.

5. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М.: Наука, 1974.

6. Дьяконов В.П., Круглов В.В. MATLAB 6.5 SP1/7/7 SP1/7 SP2 + Simulink 5/6. Инструменты искусственного интеллекта и биоинформатики. М.: СОЛОН-ПРЕСС, 2006.