CC BY
69
Вычислительные технологии
Том 18, № 3, 2013
Распознавание разрешимости интервальных
и __>к
уравнении и его приложения к анализу данных*
С. П. ШАРЫЙ, И. А. ШАРЛЯ Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: shary@ict.nsc.ru, sharia@ict.nsc.ru
Посвящается 70-летию Ю.И. Шокина
Рассматривается задача распознавания разрешимости (непустоты множества решений) интервальных систем линейных алгебраических уравнений. Для её решения развивается метод, основанный на максимизации распознающего функционала множества решений. В качестве приложения предлагается новый подход к восстановлению параметров линейных зависимостей по данным с интервальными неопределённостями ("метод максимизации согласования"). На его основе разработана методика определения выбросов в наблюдениях. Приводится сравнение метода максимума согласования с методом наименьших квадратов.
Ключевые слова: интервальная система уравнений, распознающий функционал, задача идентификации, интервальная неопределённость, метод максимума согласования.
1. Введение и обозначения
Предметом рассмотрения в работе являются интервальные системы линейных алгебраических уравнений (ИСЛАУ) вида
апх 1 + а\2Х2 + ... + а\пХп = Ьь
0,21X1 + а22Х2 + ... + а2пХп = Ь2, < . : .. . : (1)
ат1х1 + ат2х2 + . . . + атпхп Ьт
с интервальными коэффициентами а^ и интервальными свободными членами Ь, г = 1,
2, ..., т, ] = 1, 2,..., п, или, кратко,
Ах = Ь, (2)
где А = (а^) — интервальная т х п-матрица, х = (х1, х2,..., хп)т, Ь = (Ь) — интервальный т-вектор. Системы (1)-(2) понимаются как совокупности точечных систем линейных уравнений Ах = Ь с т х п-матрицами А, принадлежащими А, и т-векторами Ь из Ь.
* Работа частично финансировалась Программой Правительства России по государственной поддержке ведущих научных школ, проект № НШ-6293.2012.9.
Множеством решений интервальной линейной системы уравнений будем называть множество
Б(A, b) = { x e Rn | (3 A e A)(3 b e b)( Ax = b) }, (3)
образованное всевозможными решениями точечных систем Ax = b c A E A, b E b (см., например, [1-4]). Часто его называют также объединённым множеством решений, поскольку для интервальных уравнений и систем уравнений существуют другие множества решений, более адекватные тем или иным конкретным практическим ситуациям [3, 5-8]. В настоящей работе они не рассматриваются, и поэтому (3) называется сокращённым термином "множество решений". Если интервальная система уравнений фиксирована, то иногда для краткости множество решений будет обозначаться символом Б.
В разделе 2 работы обсуждается, как проверять пустоту или непустоту множества решений Б (A, b), как находить точку из него, а также как корректировать интервальную систему с целью достижения ею желаемых свойств. В разделе 3 рассмотрено приложение предложенной методики к задаче восстановления параметров линейных зависимостей по данным, имеющим интервальную неопределённость.
Используемая авторами система обозначений следует неформальному международному стандарту на обозначения в интервальном анализе [9]. В частности, множество всех замкнутых интервалов вещественной оси обозначено IR, а интервалы и интервальные величины представлены буквами жирного шрифта: A, B, C, ..., x, y, z. Неинтервальные (точечные) величины никак специально не выделяются. Черта снизу и сверху — X, X — означает соответственно нижний и верхний концы интервала x, так что в целом x = [x,X] = { x E R | x < x < X}. Кроме того, использованы обозначения:
mid x = 2 (x + x) — середина интервала, rad x = 1 (x — x) — радиус интервала,
|x| = max{ |x|, |x| } — абсолютное значение (модуль) интервала,
Г min{ |x|, |x|}, если 0 G x, мигнитуда интервала, (x) = < — наименьшее расстояние
0 иначе,
4 ' точек интервала до нуля.
Посредством IR помимо множества вещественных интервалов будем обозначать также классическую интервальную арифметику — алгебраическую систему, образованную интервалами вещественной оси, с операциями сложения, вычитания, умножения и деления, определёнными "по представителям", т. е. в соответствии с общим правилом
x * y = { x-ky | x e x, y e y } для * e { +, —, ■,/}.
Иными словами, результирующий интервал любой операции есть множество, образованное всевозможными результатами этой операции между элементами интервалов-операндов. Развёрнутые формулы для интервальных сложения, вычитания, умножения и деления выглядят следующим образом [1-4, 10]:
x + y = [x + У, x + y],
x — y = [x — y, x — y],
x ■ y = [min{xy,xy,xy,xy}, max{xy,xy,xy,xy}],
x/y = x ■ [1/y, 1/y] для y G 0.
Интервальный вектор определяется как вектор (столбец или строка) с интервальными компонентами. Соответственно интервальная матрица — это прямоугольная таблица, составленная из интервалов. Геометрическим образом интервальных векторов и матриц являются прямоугольные параллелепипеды с рёбрами, параллельными координатным осям пространства. Интервальные векторы часто называют также бру-сами. Для интервалов, интервальных векторов и матриц отношение включения " С" определено как включение соответствующих множеств. Операции mid и rad для интервальных векторов и матриц определяются покомпонентно и поэлементно, например, (radA)ij = rad а^. Так же покомпонентным образом применяется операция взятия модуля вектора.
Арифметические операции с интервальными векторами и матрицами являются аналогами соответствующих операций для точечного случая. В частности, сумма (разность) двух интервальных матриц одинакового размера есть интервальная матрица того же размера, образованная поэлементными суммами (разностями) операндов. Если X = (х^-) € Жтх1 и Y = (у^-) € Ж1хп, то произведение матриц X и Y есть матрица ^ = (^) € Жтхп такая, что
I
: ^ ] Хгк ук].
к= 1
Умножение интервального вектора (матрицы) на скаляр как слева, так и справа равносильно умножению всех компонент вектора (элементов матрицы) на этот скаляр.
Условимся также обозначать через т1Х топологическую внутренность множества X С Кп, т. е. множество всех точек из X, которые принадлежат множеству X вместе с некоторой своей окрестностью.
2. Теория
2.1. Разрешимость интервальных линейных систем
Говорят, что интервальная система уравнений (1)-(2)
Ах = Ь
является слабо разрешимой (слабо совместной), если существуют такие А € А, Ь € Ь, что система уравнений Ах = Ь разрешима (совместна). Интервальная система уравнений (1)-(2) называется сильно разрешимой (сильно совместной), если разрешимой является каждая система Ах = Ь с матрицей А € А и правой частью Ь € Ь [8]. Нетрудно видеть, что слабая разрешимость интервальной системы уравнений равносильна непустоте её множества решений Б (А, Ь). В настоящей работе рассматривается только слабая разрешимость интервальных систем уравнений, и поэтому будем говорить просто о разрешимости, опуская прилагательное "слабая". Отметим, что термины "сильный-слабый" в применении к свойствам интервальных уравнений и неравенств были введены в употребление И.И.Ерёминым [11, с. 93-95].
Строение множеств решений интервальных систем уравнений неплохо изучено. В частности, известно, что пересечение Б (А, Ь) с каждым ортантом пространства Кп есть выпуклое многогранное множество, для которого уравнения граничных гиперплоскостей легко выписываются по матрице и правой части ИСЛАУ (см. [3]). К примеру, на рис. 1 изображено множество решений интервальной линейной системы
Рис. 1. Множество решений интервальной системы (4)
( 4 [0, 2] [0, 2] ^
[0,2] 4 [0,2] \[0,2] [0,2] 4 у
х =
Л-1,1Л [-1,1]
На факте выпуклости и многогранности пересечения множеств решений ИСЛАУ с каждым из ортантов пространства Кга может быть основан универсальный метод распознавания разрешимости интервальных линейных систем уравнений. В самом деле, пустота или непустота пересечения Б (А, Ь) с каждым из ортантов Кга может быть выявлена путём решения некоторой системы линейных неравенств (размера 2т х п), например, с помощью хорошо разработанных методов линейного программирования. В целом же распознавание этим способом множества решений ИСЛАУ и нахождение его точки потребуют не более чем 2П (количество ортантов пространства Кга) решений систем линейных неравенств. Ясно, что практически этот рецепт работоспособен лишь при малом числе п неизвестных задачи.
Трудоёмкость описанного подхода не может быть существенно улучшена, что следует из принципиальной труднорешаемости рассматриваемой задачи, т. е. того, что она КР-трудна [8, 12, 13]. Таким образом, в общей ситуации исследование пустоты/непустоты множества решений ИСЛАУ и нахождение точки из него является весьма непростым делом.
2.2. Метод распознающего функционала
Для распознавания разрешимости интервальных систем линейных алгебраических уравнений в работе развивается техника распознающих функционалов, разработанная ранее в [14-17].
Введём в рассмотрение функционал Uss : х IRmxn х IRm ^ R, определяемый
как
Uss (x, A, b) = min < rad bi + > (rad a^) |xj | —
1 <T iCm I ' ^
j=i
mid bi — ^^ (mid aij) j=i
ij) xj
где х — вещественный п-вектор, А — интервальная т х п-матрица, Ь — интервальный т-вектор. Там, где подразумевается, что интервальные переменные фиксированы и необходимо подчеркнуть, что функционал Шя рассматривается только в зависимости от переменной х, будем писать Ивв (х), опуская упоминание интервальной матрицы А и интервального вектора Ь.
Теорема 1. Пусть А — интервальная тхп-матрица, Ь — интервальный т-век-тор. Тогда принадлежность точки х € Ега множеству решений Е(А, Ь) интервальной линейной системы Ах = Ь равносильна неотрицательности в х функционала иве:
x е ^ (A, b)
Uss (x, A, b) > 0.
Иными словами, множество решений Е(А, Ь) интервальной системы линейных уравнений Ах = Ь является лебеговым множеством { х € Ега | Ивв (х, А, Ь) > 0 } функционала Ивв по переменной х.
Доказательство. Напомним, что для любого вектора х € Кга интервальное мат-рично-векторное произведение Ах совпадает с множеством произведений "по представителям" { Ах | А € А }. По этой причине непустота пересечения Ах и Ь означает существование таких А € А, Ь € Ь, что Ах = Ь. Это в свою очередь эквивалентно х € Е(А, Ь). В целом
х € Е(А, Ь) ^^ Ах П Ь = 0. (6)
Данная равносильность известна в интервальном анализе как характеризация Бека для множества решений интервальных линейных систем уравнений [3, 4, 18].
Далее, два интервала р и д вещественной оси имеют непустое пересечение тогда и только тогда, когда
р — midд| < rаdр + radд.
Поэтому непустота пересечения двух интервальных брусов Ах и Ь имеет место тогда и только тогда, когда
mid(Ax)i — mid bi| < rad(Ax)i + rad bi
i = 1, 2,
, m.
Коль скоро mid(Ax) = (mid A) x и rad(Ax) = (rad A) |x| (см. [1, 3, 4]), выписанному неравенству можно придать следующий вид:
| ((midA) x)i — mid bi| < ((radA) |x|)i + rad bi, i = 1, 2, или в векторной форме
| (midA) x — mid b| < (radA) |x| + rad b.
, m,
(7)
Это соотношение называется неравенством Оеттли —Прагера [3, 4, 8].
Перепишем неравенства (7) в равносильном виде
rad 5» + ((radA) |х|)» — I mid 5» — ((midA) > 0, i = 1, 2,..., m,
или развёрнуто
rad 5» + ^^ (rad a^-) | — j=i
mid 5» — ^^ (mid a^-) j=i
j)x j
> 0, i = 1, 2,..., m.
Далее посредством взятия минимума по г можно свернуть т штук условий (8) в одно, заключая, что точка х принадлежит множеству решений Б (А, Ь) в том и лишь в том случае, если
mm < rad 5» +> (rad a„) \xA— l<»<m I
j=i
mid 5» — ^^ (mid a»,) j=i
ij)x j
0.
Это и требовалось доказать. □
Итак, функционал Шя знаком своих значений позволяет "распознавать" принадлежность точки множеству Б (А, Ь). По этой причине будем называть его распознающим функционалом множества решений Б (А, Ь). Выражения (8), т.е.
^i(x) = rad 5» + ^^ (rad a^) | — j=l
mid 5» — ^^ (mid aj) j=l
j)x j
, i =1, 2,..., m,
минимумом или нижней огибающей которых является функционал Шя, будем называть его образующими.
Ясно, что доказанная теорема останется в силе, если вместо Шя взять функционал
Шв 7(х, А, Ь) =
mm < y» I rad 5» + ^^ (rad a^) |
1<»<m j=i
mid 5» — ^^ (mid a j) j=i
»j)x j
в котором образующие имеют какие-то положительные веса 7^, г = 1, 2,..., т. Далее в разделах 2.3 и 3.3 будет показано, что в некоторых ситуациях более уместно применение именно такой модифицированной формы распознающего функционала. Другая интересная форма распознающего функционала Шя имеет вид
Uss (я, A, 5) = min < rad(Ax — 5)» — 1 mid(Ax — 5)J 1<»<т I J
:iq)
Эту форму И. Рон сообщил авторам настоящей статьи в процессе её подготовки. Приведём вывод представления (10) из характеризации Бека.
Непустота пересечения Ах П Ь = 0 равносильна включению 0 Е Ах — Ь. Но интервальный брус содержит нулевой вектор тогда и только тогда, когда каждая его компонента содержит нуль. Интервал содержит нуль тогда и только тогда, когда его радиус не меньше модуля середины. Как следствие, получаем (10).
То, что (10) и (5) задают один и тот же функционал, следует из равенств
rad(Ax — b) = rad(Ax) + rad b = rad b + (radA) |x|, | mid(Ax — b) | = | mid(Ax) — mid b | = | mid b — (midA) x |.
Функционал Uss (x, A, b) очевидно непрерывен по любому своему аргументу. Из непрерывности по x вытекает общеизвестный результат: множество решений интервальной системы линейных алгебраических уравнений замкнуто.
Дальнейшие свойства функционала Uss сформулируем в виде последовательности предложений.
Предложение 1. Функционал Uss (x, A, b) является непрерывным по Липшицу относительно любого из своих аргументов.
Это следует из того, что липшицевыми являются все операции и функции, из которых составлено выражение для функционала Uss.
Предложение 2. Функционал Uss — вогнутый по переменной x в каждом ор-танте пространства Rn. Если в матрице A интервальными являются l столбцов с номерами из множества J = (ji, j2,..., ji}, l < n, а остальные столбцы матрицы — целиком точечные (неинтервальные), то Uss (x, A, b) вогнут по x на объединениях нескольких ортантов, точнее, на каждом из 21 множеств вида
{ x е | xj ^ 0,j е J },
где означает какое-либо из отношений "< или ">".
Доказательство. Достаточно провести его для всех образующих
^i(x) = rad bi + ^^ (rad a^) |xj | j=i
mid bi — ^^ (mid aij-) j=i
ij) x j
i = 1, 2,..., m,
поскольку Uss является их нижней огибающей.
В пределах одного ортанта из Кп, когда знаки компонент х^- вектора х остаются постоянными, функции
п
rаd Ь + ^^ (rad а^) 1 (11)
¿=1
являются линейными по х. Кроме того, функции
mid bi — ^^ (mid aij) j=i
ij) xj
глобально вогнуты. Как следствие, в пределах одного ортанта все ^¿(х), г = 1, 2,... , т, будучи суммой линейной и вогнутой функций, также вогнуты относительно х.
В общем случае провести выполненные выше рассуждения нельзя, так как выраже-
п
ние £ (rаd а^^) 1 не является линейной функцией от х за пределами одного ортанта ¿=1
пространства Кп. Но если для некоторого индекса к Е {1, 2,..., п} все а^^ — точечные (неинтервальные), то rad а^^ = 0, г =1, 2,... , т, и выражения (11) линейны для любых значений .
Итак, если в А точечным является целый к-й столбец, так что radа^ = 0, г = 1, 2,... , т, то и все функции ^¿(х) (а вместе с ними и Шя) будут вогнутыми на множествах { х € Ега | х1 ^ 0,... , хк-1 ^ 0, хк+1 ^ 0,... , хп ^ 0 }, каждое из которых является объединением двух ортантов в Кга. Обобщение этого рассуждения на случай нескольких точечных столбцов в А очевидно. □
Предложение 3. Функционал Шв — кусочно-линейный, т. е. его график составлен из конечного числа кусков гиперплоскостей.
Доказательство. Так как |а| = шах{а, —а}, то в пределах каждого отдельного ортанта пространства Кга, где знаки компонент х постоянны, имеем
^j(x) = rad b + ^^ (rad aj) sgn Xj ■ X^ -
j=i
— max < mid b — ^^ (mid aj) Xj, -mid b + ^^ (mid aj) Xj
I j=i j=i
i = 1, 2,... , m. Поэтому в любом ортанте
Uss (x) = min ^i(x) =
1< i<m
rad bi + ^^ (rad aij) sgn Xj ■ Xj — mid bi + ^^ (mid aij) j=i j=i
rad bi + ^^ (rad aij) sgn Xj ■ Xj + mid bi — ^^ (mid aij) Xj j=i j=i
in
— bi + > (rad aij sgn Xj + mid a^) Xj, z—/ v /
j=i
n
bi + ^^ (rad aij sgn Xj — mid a^) Xj j=i
= min min ( — b + a-X, bi — afX) , (12)
i<i<m - ' - - J '
где a- и a- — какие-то вершины интервальной вектор-строки (aii, ai2,..., ain), определяемые как
(a-) j = rad aij sgn Xj + mid aij,
(af) . = rad aij sgn Xj — mid aij, j = 1, 2,..., n.
'з
Выражения в фигурных скобках из (12), по которым берутся минимумы, задают линейные функции. Таким образом, функционал Шв (х) в каждом ортанте представляет собой минимум конечного числа линейных функций и потому в целом кусочно линеен. □
В общем случае функционал Шв (х, А, Ь) может быть неограниченным сверху по переменной х на всём пространстве Кга, как показывает следующий пример. Для одномерного интервального уравнения [—1, 3] х = [1, 3] множеством решений является [—то, —1] и [1/3, Распознающий функционал задаётся выражением
Шв (х) = 1 + 2 |х| — |х — 2| (13)
и очевидно неограничен. Его график представлен на рис. 2.
Рис. 2. График неограниченного распознающего функционала (13)
Предложение 4. Если множество решений Е(А, Ь) ограничено, то функционал Шв (х) ограничен сверху.
Доказательство. Если множество решений Е(А, Ь) пусто, то в силу теоремы 1 функционал Шб (х) ограничен сверху нулём.
Если множество решений ИСЛАУ непусто и ограничено, то распознающий функционал Шв (х), будучи непрерывной функцией, является ограниченным на множестве решений и, как следствие, ограниченным сверху на всём □
В целом выявление необходимых и достаточных условий глобальной ограниченности для Шя сверху является интересной открытой задачей. Отметим, что если распознающий функционал неограничен сверху, то множество решений соответствующей интервальной системы уравнений непусто.
Предложение 5. Если функционал Шб (х) ограничен сверху, то он достигает конечного глобального максимума на всём пространстве Ега.
Доказательство. Достаточно показать, что распознающий функционал Шя достигает конечного максимума в каждом ортанте О С
Далее понадобится понятие подграфика функционала Шя. Это множество, определяемое как
hypUss = {(х,г) е х к | г < Шб (х)}.
Если О — ортант в то в силу Предложений 2 и 3 пересечение подграфика с множеством Ох К пространства Кга+1 является выпуклым многогранным. Как всякое выпуклое многогранное множество, оно представимо (см. [19]) в виде комбинации конечного числа точек (^ , ^), ^ е ^ е К, ] = 1, 2,...,р, и направлений (^, в,), ф е б^- е К, ] = 1, 2,... , д, так что
p q
hypUss П (Ox R)=^ Aj (Sj ,tj) + J] ^ (dj Aj > 0,^- > 0,J]Aj = 1 |> . (14
j=i j=i
p
j=i
Из ограниченности функционала Шя сверху следует, что в представлении (14) должно быть е^ < 0, в противном случае в множестве hypUss П (О х К) окажутся точки со сколь угодно большой (п + 1)-й координатой. Как следствие,
max Uss (x) = maxi (n + 1)-я координата точек из hyp Uss П (Ox R) }
x€O
p q p
Y1 Aj^ + Y1 ^ej Aj >0, ^ >Aj =
^frax \ / Aj tj
p
max ( ^^ Aj tj
Aj > 0, > Aj = 1 > = max tj,
i<j<p
j=i
j=i
и этот максимум достигается в той точке (sj,tj) из комбинации (14), которая имеет наибольшую (n + 1)-ю координату. □
Отметим, что это доказательство почти дословно совпадает с доказательством Предложения 4 из [17] о том, что другой распознающий функционал Uni (x) всегда достигает глобального максимума.
Предложение 6. Если Uss (x, A, b) > 0, то x — точка топологической внутренности int Е(A, b) множества решений, т. е. принадлежит множеству Е(A, b) вместе с некоторой своей окрестностью.
Доказательство. Поскольку функционал Uss (x) непрерывен, то множество Y = {y G | Uss (y) > 0 } — открытое. Кроме того, по условию x G Y, а из теоремы 1 следует, что Y Ç Е(A, b). По этой причине int Е(A, b) = 0 и x G int Е(A, b). □
Доказанное свойство распознающего функционала позволяет использовать его для исследования того, принадлежит ли точка внутренности множества решений. Это может иметь особую важность при нахождении телесной внутренней оценки множества решений вокруг точки-центра по методике, описанной в [3, 20].
Предложение 7. Пусть интервальная линейная система уравнений Ax = b такова, что её расширенная матрица (A, b) не содержит строк, все элементы которых имеют нулевые концы. Тогда из принадлежности x G int( Е(A, b) П О), где О — ор-тант пространства Rn, следует строгое неравенство Uss (x, A, b) > 0.
Доказательство. Оно показывает, что если Uss (x) = 0, то равенство нулю функционала Uss имеет место также для всех x G Е ПО. Но это невозможно при наложенных на систему условиях.
Пусть max Uss (x) достигается в точке u G Е ПО. Из выпуклости множества Е ПО и
условия x G int (ЕПО) следует, что точка x является внутренней для некоторого отрезка прямой [u,y], который соединяет u и y и лежит в Е П О, в силу чего x = Au + (1 — A)y, 0 < A < 1. Поскольку функционал Uss вогнут на О, то
Uss (x) > A Uss (u) + (1 — A) Uss (y), (15)
где A > 0, (1 — A) > 0. Если предположить, что Uss (x) = 0, то, принимая во внимание неравенства Uss (u) > 0 и Uss (y) > 0, можно заключить, что выражение (15) может быть справедливым лишь при Uss (u) = Uss (y) = 0. Но Uss (u) = mmaxUss, и тогда
функционал Uss должен быть тождественно нулевым всюду на Е П О.
При доказательстве Предложения 3 было показано, что в каждом ортанте функционал Uss (x) является минимумом конечного числа линейных функций:
Uss (x) = min min { — b + aix, b — afx} ,
1<i<m
где a', ai' — какие-то вершины интервальной вектор-строки (ai1, ai2,..., ain). Зануление функционала Uss (x) на телесном множестве Е (A, b) ПО означает, что хотя бы одна из линейных функций вида (—b+aix), (b — ai'x), i = 1, 2,..., m, должна быть нулевой. Это противоречит тому, что расширенная матрица (A, b) не содержит строк, все элементы которых имеют нулевые концы.
Полученное противоречие приводит к заключению, что должно иметь место строгое неравенство Uss (x) > 0. □
Основываясь на полученных выше результатах, естественно приходим к конструктивному способу исследования разрешимости интервальных линейных систем уравнений, который мы называем методом распознающего функционала. Метод заключается в том, что для интервальной системы уравнений Ax = b решается задача безусловной максимизации по x распознающего функционала Uss (x, A, b) и в зависимости от результата делается вывод о разрешимости ИСЛАУ и свойствах её множества решений. Рассмотрим сначала случай, когда удалось найти глобальный максимум
M = max Uss (x, A, b), xei"
который достигается в точке т Е ln. Тогда
— если M > 0, то т Е £(A, b) = 0, т.е. интервальная линейная система Ax = b разрешима и т лежит в множестве решений;
— если M > 0, то т Е int £(A, b) = 0, и принадлежность точки т множеству решений устойчива к малым возмущениям A и b;
— если M < 0, то £ (A, b) = 0, т.е. интервальная линейная система Ax = b неразрешима.
Кроме того, знание конкретной величины конечного максимума M = max Uss (x, A, b)
позволяет более тонко исследовать интервальную систему уравнений.
Рассмотрим теперь случай, когда не удалось найти глобальный максимум распознающего функционала Uss (x). Тогда информация, полученная в процессе его поиска, всё-таки может быть полезной при исследовании разрешимости системы Ax = b:
— если установлена неограниченность Uss (x) сверху, то интервальная система уравнений Ax = b разрешима;
— если найдена такая точка т, что Uss (т) > 0, то т Е £(A, b) = 0, т. е. интервальная линейная система Ax = b разрешима и т лежит в множестве решений;
— если найдена такая точка т, что Uss (т) > 0, то т Е int £(A, b) = 0 и принадлежность точки т множеству решений устойчива к малым возмущениям A и b.
При отсутствии перечисленных данных вопрос о разрешимости интервальной системы линейных уравнений Ax = b остаётся открытым.
2.3. Коррекция интервальной системы уравнений
Специфический вид выражения для распознающего функционала Uss даёт возможность легко прогнозировать, как будут изменяться его значения при изменении интервалов в матрице A и правой части b интервальной системы линейных алгебраических уравнений Ax = b.
Заметим, что величины rad b^ входят слагаемыми во все образующие ^¿(x), нижней огибающей которых является распознающий функционал Uss. Поэтому одновременное изменение всех rad bi на одну и ту же величину приводит к тому, что значение распознающего функционала изменится на ту же величину. В частности, если C > 0
и e = ([-1,1],..., [-1,1])т, то для системы Ax = b + Ce с уширенной на Ce правой частью имеем Uss (x, A, b + Ce) = Uss (x, A, b) + C и, следовательно,
max Uss (x, A, b + Ce) = max Uss (x, A, b) + C.
же!п xei"
Если множество решений ^ (A, b) пусто, так что
M = max Uss (x, A, b) < 0,
то при увеличении радиуса всех компонент правой части на C > |M | это множество сделается непустым, поскольку max Uss (x, A, b + Ce) в силу (16) станет уже неотрицательным.
Рассмотрим теперь случай, когда требуется сделать уширение компонент правой части не одинаковым, а пропорциональным положительным весам xi > 0, i =1, 2,... , m. Тогда следует найти MY = max Uss Y(x, A, b) для модифицированного распознающего
функционала (9) с Yi = к—1, i =1, 2,... , m, и увеличить радиус каждого b на XiC, где C > M|.
Возможность коррекции интервальной линейной системы путём изменения её матрицы основана на следующих соображениях. Увеличение радиуса всех элементов матрицы A = (aj ) на положительную величину C равносильно прибавлению к A интервальной матрицы CE, где E — m х n-матрица, все элементы которой суть [-1,1]. Если в точке x = 0 распознающий функционал принимает значение U = Uss (x, A, b) < 0, то для интервальной системы с уширенной матрицей (A + C E), которая образована элементами (aj + C [-1,1]), имеем
Uss (x, A + CE, b) =
min -j rad b + ^^ (rad a j + C) |xj |
|xj I -
j=i
mid b — ^^ (mid a j)
j) x j
j=i
min ^ rad b + ^^ (rad a j) |xj | + ^^ C |
|xj | -
j=i
j=i
mid b — ^^ (mid a j )
-ij ) xj
j=i
min < rad bi + ^^ (rad aij ) |xj | j=i
mid bi — ^^ (mid aij ). j=i
+ C £ |xj |
j=i
коль скоро выражения С ^ 1 не зависят от индекса г. Поэтому если увеличение радиусов всех элементов исходной матрицы А на величину С таково, что
C £ |xj| > |U|,
j=i
то получим Uss (х, А+СЕ, Ь) > 0. Следовательно, точка х будет принадлежать непустому множеству решений Б (А+С Е, Ь). Как и для правой части, введение положительной матрицы весовых множителей позволяет добиться корректировки исходной системы путём неодинакового уширения элементов матрицы А.
2.4. Сравнение распознающих функционалов
Ранее в работах [14, 17] был введён другой распознающий функционал множества решений интервальной системы линейных алгебраических уравнений. Он называется Uni и задаётся выражением
Uni(x, A, b) = min
1< i<m
rad b — ^ mid b — ^^ aj xj ^ |
'17)
Принадлежность точки x G множеству решений S (A, b) интервальной системы линейных уравнений Ax = b равносильна неотрицательности в x функционала Uni. Свойства Uni близки к свойствам Uss, но в целом эти функционалы заметно отличаются друг от друга. Иногда использование Uss является более выгодным, чем функционала Uni, а иногда наоборот.
Распознающий функционал Uss в том виде, как он даётся выражением (5), несколько проще, чем Uni, так как для вычисления Uss вообще не используются интервальные арифметические операции, достаточна только обычная вещественная арифметика. В выражении для Uss нет такой относительно сложной операции, как взятие мигнитуды. С другой стороны, функционал Uni всегда ограничен сверху и достигает конечного максимума, чего в отношении Uss в общем случае сказать нельзя.
Исследуем, как соотносятся между собою значения функционалов Uss и Uni. Для любого i = 1, 2,... , m имеют место соотношения
mid b — ^^ a
j=i
ij xj
mid bi — ^^ (mid aij + [—1,1] rad a j=1
ij xj
mid bi — ^^ (mid aij-) xj — ^^ [—1,1] (rad aij- )
j=1 j=1
ij )x j
в силу дистрибутивности умножения по сложению для точечных общих множителей х^ (см. [1-4])
> / mid bi — ^^ (mid aij ) ;
j=1
^ 1] (rad aij );
j=1
в силу свойства мигнитуды (p + g) > (p) — |q| (см. [3, 4])
mid bi — ^^ (mid aij) j=i
ij) x j
J^(rad aij) |xj1 j=i
так как все [-1,1] (rad а^) х^ — уравновешенные (симметричные относительно нуля интервалы). Поэтому
rad bi — ( mid bi — ^^ aij xj ) < rad bi + ^^ (rad aij ) |xj |
j=i
j=i
mid bi — ^^ (mid aij) ; j=i
Коль скоро это верно для всех i = 1, 2,..., m, то
Uni (x) < Uss (x)
при любых x G Rn, т. е. значения функционала Uss всегда не меньше значений функционала Uni.
Если правая часть ИСЛАУ имеет хотя бы одну точечную компоненту bi с rad bi = 0, то, учитывая неравенство
/ n \
/ mid bi — ^^ aij xj ) > 0,
можно видеть, что функционал Uni при непустом и телесном (с непустой внутренностью) множестве решений в качестве максимума своих значений может иметь плато нулевого уровня, тогда как у функционала Uss максимум достигается при этом в "острой" вершине графика, строго большей нуля. Последнее обстоятельство более благоприятно при численном нахождении этого максимума, а также при распознавании разрешимости интервальных линейных систем, когда максимум распознающего функционала необходимо сравнивать с нулём.
Для иллюстрации на рис. 3 изображены графики функционалов Uni и Uss для интервальной линейной системы
( ЙИ) x = ( |-0'. (18)
В целом, по-видимому, нельзя однозначно утверждать, что функционал Uss лучше или хуже распознающего функционала Uni. Эти функционалы просто разные, ведут себя в различных ситуациях по-разному и, как следствие, предназначены для различных целей. В частности, функционал Uni "сильнее" учитывает правую часть ИСЛАУ и потому может быть более полезным, к примеру, в тех прикладных задачах, где правая часть (соответствующая неопределённостям на выходе в задачах идентификации, которые рассмотрены в следующем разделе) играет решающую роль.
2.5. Обобщения на случай нелинейных систем уравнений
Рассмотрим общий случай нелинейной интервальной системы уравнений вида
F (a, x) = b,
где F = (Fi(a,x),F2(a,x),...,Fm(a,x))T, x = (xi,x2,...,xn)T, Fi : R1 х Rn ^ R — заданные функции, a = (ai, a2,..., a^, b = ( bi, b2,..., bm) — интервальные векторы. Множеством решений этой системы называется множество
Б (F, a, b) = { x G Rn | (3 a G a)(3 b G b)( F (a,x) = b ) }, (19)
образованное всевозможными решениями систем F (a, x) = b для a G a, b G b. В отношении этого множества естественно возникает задача исследования его непустоты. Будем далее предполагать, что все Fi являются непрерывными функциями своих аргументов.
Аналогично линейному случаю можно попытаться ввести распознающий функционал множества решений, который мог бы выделять точки множества решений (19) с помощью знака своих значений. Тогда задача исследования множества решений на непустоту сведётся к максимизации этого распознающего функционала. Но, конструируя этот функционал в общем нелинейном случае, мы сталкивается с большими трудностями при определении области значений многомерного отображения F : R1 х Rn ^ Rm, так как при произвольных F эта область значений может иметь сложную форму и не быть брусом в Rm. Если обозначить как F (a, x) область значений отображения F (a,x)
Рис. 3. Графики распознающих функционалов Uni и Uss для системы (18)
по a G a, то нелинейным аналогом характеризации Бека (6) может служить эквивалентность
x G S (F, a, b) ^^ F (a, x) П b = 0.
Но в общем случае почти невозможно придать ей какую-либо удобную аналитическую форму.
Тем не менее, если наложить на отображение F(a,x) условие из [5, 7] что каждый параметр aj, j = 1, 2,... , /, входит не более чем в одну из компонент Fi, i = 1, 2,... , m,
то область значений F (a, x) является прямым произведением областей значений отдельных компонент Fj(a,x), т.е.
F(a,x) = Fi(a,x) х F2(a,x) х ■ ■ ■ х Fm(a,x),
где Fj(a,x) = { Fj(a,x) | a G a } — области значений функций Fj по переменной a G a. Соответственно, F (a, x) П b = 0 равносильно тому, что
Fj(a, x) П b = 0 для всех i = 1, 2,..., m.
Тогда можно ввести распознающие функционалы для множества решений (19) —
UniF(x, a, b) = imip { rad b — ^ mid b — Fj(a, x) ^ |
и
UssF(x, a, b) = min
1<i<m
свойства которых аналогичны свойствам функционалов Uni и Uss для линейного случая. Ясно, что полезность этих функционалов зависит от функций Fj(a, x), i = 1, 2,..., m, и эффективности вычисления их областей значений. В самом общем случае эта задача непроста. Но если выражение для Fj(a,x) имеет по одному вхождению каждой переменной, то её область значений вычисляется относительно просто, как естественное интервальное расширение [2, 3, 10].
i rad b + rad Fj(a, x) — mid b — mid Fj(a, x)
3. Идентификация в условиях интервальной неопределённости
В качестве примера практического применения развитой выше методики исследования разрешимости интервальных систем уравнений рассмотрим задачу идентификации параметров системы в условиях интервальной неопределённости.
3.1. Задача идентификации
Разберём сначала неинтервальный случай. Пусть имеется статический (не зависящий от времени) объект, входы которого описывается вектором a = (a1, a2,... , an) G Rra, на выходе имеем величину b G R (рис. 4), а зависимость вход-выход является линейной, имеющей вид
b = ai xi + a2x2 + ... + araxra (20)
с некоторыми постоянными вещественными коэффициентами Xi, Х2, . . . , Xn. Задача идентификации (называемая также задачей восстановления зависимостей) — это задача
al a2
Объект
b
a
n
Рис. 4. Структурная схема объекта идентификации
определения (оценивания) значений на основе данных о входе и выходе объекта, т. е. по ряду соответствующих друг другу значений a и b.
Каждое наблюдение (измерение) входов и выходов объекта порождает соотношение, которое связывает искомые . Если серия наблюдений входа-выхода является "достаточно представительной", то можно попытаться решить получающуюся систему уравнений относительно неизвестных ж^ и найти их значения. Восстанавливаемую функциональную зависимость (20) обычно называют регрессией величины b по величинам ai, a2, ..., an, а её график — регрессионной линией b по a1, a2, ..., an.
Для удобства условимся обозначать результат i-го измерения входов такого объекта через (ai1, ai2,..., ain), а выхода — через b¿, i = 1, 2,..., m. Тогда решение задачи идентификации — это вектор ж = (ж1, ж2,..., жп)т, который удовлетворяет всем соотношениям вида (20), получающимся в результате отдельных наблюдений. Иными словами, результат идентификации ж = (ж1,ж2,... , жп)т является решением системы уравнений
ai^i + a^2 +... + a^n = b1,
a2^1 + a22Ж2 + ... + a2raЖ„ = b2,
. . .. . : (21)
am1Ж1 + am2Ж2 + ... + amn жп bmo
где m — общее количество наблюдений (измерений).
На практике эта ясная и стройная схема почти не работает, поскольку система уравнений (21), как правило, решений не имеет. Это вызвано тем, что данные измерений (наблюдений) содержат случайные ошибки, которые делают невозможными точные равенства (20). Дополнительно неразрешимости системы (21) способствует то обстоятельство, что она обычно переопределена: измерений стремятся получить как можно больше и заведомо не меньше количества неизвестных параметров, подлежащих определению. Кроме того, назначая избранный вид зависимости между входными и выходными переменными, мы можем также совершать идеализацию, желая, например, выявить только "главную часть" функциональной зависимости, которая связывает b и a1, a2, ..., an. Как следствие, для системы (21) необходимо искать решение в каком-то обобщённом смысле. Чаще всего в качестве такого решения ищут точку (или точки), минимизирующую норму невязки, т. е. разности между правой и левой частями системы (21).
3.2. Идентификация по данным с интервальной неопределённостью
Рассмотрим более сложную практическую ситуацию, когда входы и/или выходы объекта не известны точно и доступны лишь их интервальные оценки
a¿j е a¿j = [ , a^] и b е b =[ 5,, b,]. (22)
Эта неопределённость может быть следствием ошибок измерений или вытекать из принципиальных трудностей в точном определении значений некоторых величин, а применение других моделей неопределённости иногда оказывается просто невозможным (см., к примеру, подробные комментарии к постановке задачи в [22-27]). В этих новых условиях приходим к необходимости рассматривать задачу идентификации при интервальных данных (22). Впервые подобную задачу поставил Л.В. Канторович [28] (он предполагал наличие неопределённости лишь в выходных переменных).
Единого универсального понимания решения подобной интервальной идентификационной задачи быть не может, но общая естественная идея — это требование, по возможности, наиболее полного "согласования", в том или ином смысле, решения идентификационной задачи с экспериментальной информацией. Понятно, что большее согласование с данными имеет тот набор параметров, которому соответствует регрессионная линия, имеющая "меньшее отклонение" от экспериментальных данных. При этом абсолютное согласование достигается, очевидно, в случае точного прохождения регрессионной линии через все точки измерений.
Сказанное выше мотивирует следующее широко принятое определение (которое иногда подразумевается неявно): набор параметров XI, х2, ... , жга объекта, описываемого (20), называется согласованным (совместным) с интервальными экспериментальными данными (ал, аг2, ..., агп), Ьг, г = 1, 2,... ,т, если для каждого г (т.е. для каждого наблюдения) в пределах измеренных интервалов существуют точечные представители ац € ал, аг2 € аг2, ..., агп € агп и 6г € Ьг такие, что выполнено отношение
«¿1X1 + а^2Ж2 + ... + аг„ж„ = 6г. (23)
Если А = (а^) — интервальная т х п-матрица, составленная из т результатов измерений входов, Ь = (Ь1, Ь2,..., Ьт)т — интервальный вектор т измерений выходов, то семейство всех векторов параметров, согласующихся с интервально заданными экспериментальными данными, может быть представлено в виде
{ х € Ега | (ЗА € А)(36 € Ь)(Ах = 6)},
т.е. как множество решений всевозможных точечных систем Ах = 6 с А € А, 6 € Ь. Специалистами по теории идентификации это множество часто называется информационным множеством [29, 30], (апостериорным) множеством возможных значений параметров, множеством допустимых значений параметров [24] и т. п. Это множество является в точности множеством решений Б (А, Ь) интервальной линейной системы уравнений Ах = Ь, как оно было определено выше в (3).
Геометрически поиск параметров зависимости означает выбор гиперплоскости, задаваемой значениями х1, х2, ..., хп, которая проходит через все брусы (аг1, аг2, ..., агп, Ьг), г = 1, 2,..., т, соответствующие результатам измерений. В двумерном случае наглядной иллюстрацией этому может служить рис. 5.
Рис. 5. Задача выбора параметров, согласующихся с данными, равносильна проведению гиперплоскости через брусы неопределённости данных
Несмотря на простоту и естественность выписанного определения результата идентификации, оно оказывается весьма жёстким и способно приводить к парадоксу. Впервые на него указал Е.З. Демиденко в [31] (см. также [32]), и его суть может быть охарактеризована фразой "чем лучше, тем хуже".
Дело в том, что присутствие неопределённости в данных — нежелательное явление, которое искажает истинную картину реальности. Уменьшение этой неопределённости, т. е. сужение интервалов данных, должно рассматриваться как положительный фактор, так как оно способствует большей точности исследуемых моделей и получаемых на их основе результатов. С другой стороны, при более широких интервалах исходных данных множество решений интервальной системы также является более широким, и поэтому появляется больше возможностей для выбора из него параметров модели, чем в случае узких интервальных данных, когда множество решений может стать вообще пустым. Итак, чем выше точность исходных данных, чем меньше их интервальная неопределённость, тем хуже для оценивания параметров. И наоборот, чем шире интервальная неопределённость, чем меньше мы знаем о точных значениях интересующих нас величин, тем лучше для процесса оценивания параметров и тем более богатый набор результатов можно получить.
Причина парадокса Е.З. Демиденко заключается в том, что результат идентификации опирается на непустоту информационного множества задачи, изменяющегося в зависимости от данных таким образом, который при ближайшем рассмотрении кажется парадоксальным. Для лучшего понимания природы возникающего затруднения следует отметить, что в предельном случае, который охватывается традиционной статистикой и оперирует с точечными (неинтервальными) данными, абсолютное и полное согласование с ними искомых параметров является нетипичным и разрушается при сколь угодно малых возмущениях в этих данных (более точно, мера множества таких ситуаций равна нулю). Само понятие информационного множества теряет при этом содержательный смысл, так как в этих условиях оно по большей части пусто. Но принципиальная особенность идентификации в случае интервальных данных состоит в том, что множество абсолютно согласованных параметров часто может быть непустым и даже телесным (иметь непустую внутренность). При этом оно устойчиво к возмущениям в данных.
В целом подходы, которые являются интервальными расширениями традиционных точечных методов (метода наименьших квадратов, наименьших модулей, чебышёвского выравнивания и др.), неудовлетворительно обрабатывают ситуации, когда информационное множество непусто. Поскольку для любой точки этого множества невязка точечных систем уравнений нулевая, то множество аргументов минимума целевой функции — целое плато нулевого уровня. Это плохо, в частности, тем, что исследование вопроса о принадлежности той или иной точки информационному множеству задачи требует сравнения с нулём. В условиях приближённого характера вычислений на ЭВМ подобная операция некорректна.
С другой стороны, подходы, ориентированные на работу с непустым информационным множеством, сталкиваются с затруднениями, когда информационное множество пусто. В подобных случаях идентификация, как правило, требует какого-то дополнительного приёма для определения параметров. Так происходит, в частности, в методе центра неопределённости, развиваемом Н.М. Оскорбиным, С.И. Жилиным и др. [24, 33-36].
Резюмируя краткий обзор методов идентификации в условиях интервальной неопределённости, можно сказать, что необходим подход, который бы единообразно обраба-
тывал ситуации, возникающие как для пустого, так и непустого информационного множества задачи.
3.3. Метод максимума согласования
Общий подход к идентификации параметров состоит в том, чтобы выбрать количественную меру "согласования" данных с параметрами объекта, и затем искомой оценкой можно будет, к примеру, взять ту точку пространства параметров, в которой эта согласованность максимальна.
Какой именно можно взять количественную характеристику согласования данных и параметров модели? Существуют естественные требования, которым "мера согласования" должна удовлетворять. При непустом информационном множестве она должна быть неотрицательной для тех точек из этого множества, на которых согласование в самом деле достигается. При этом для точек внутренности множества решений мера согласования должна быть очевидно не меньшей (или даже большей), чем на границе. Отсутствие согласования (противоречивость данных) можно считать "отрицательным согласованием", поэтому в целом мера согласования может быть каким-либо функционалом, принимающим вещественные значения из R.
Свойства распознающего функционала Uss, исследованные в разделе 2.2, неплохо удовлетворяют сформулированным выше условиям, так что Uss вполне подходит на роль искомой меры согласования. При этом отрицательные значения Uss можно понимать как "меру несогласования". Взвешенная форма распознающего функционала UssY (9) способна учитывать различную ценность рассматриваемых наблюдений. Как следствие, приходим к методу оценивания параметров, который можно назвать "методом максимизации согласования" (или методом максимума согласования): искомой оценкой параметров берётся точка arg max Uss , в которой достигается наибольшее значение распознающего функционала Uss информационного множества. При этом
— если max Uss > 0, то эта точка лежит в непустом множестве параметров, согласующихся с данными;
— если max Uss < 0, то множество параметров, согласующихся с данными, пусто, но в этой точке минимизируется "несогласованность" с данными.
Из полученных в разделе 2.3 результатов следуют дальнейшие содержательные интерпретации метода максимума согласования. Так, если множество параметров, согласующихся с данными, пусто, то arg max Uss — первая точка, которая появится в этом множестве при равномерном уширении вектора правой части относительно его середины. Таким образом, согласно методу максимизации согласования при пустом информационном множестве выбирается в качестве значений искомых параметров точка, на которой минимизируется увеличение интервальной неопределённости в выходных переменных, делающее это множество непустым. Подобный принцип согласования противоречивых данных был положен в основу метода центра неопределённости [24, 33-35]. В методе максимизации согласования он автоматически получается как следствие из общего правила выбора оцениваемых параметров.
Другой подход к интерпретации результатов метода максимума согласования состоит в том, что arg max Uss соответствует регрессионной линии, которую необходимо минимально уширить до такой "регрессионной полосы", что она уже будет иметь согласование с данными, т.е. пересекать все брусы неопределённости данных (рис. 6).
Рис. 6. "Регрессионная полоса" вместо регрессионной линии
Отметим также, что метод максимизации согласования свободен от рассмотренного выше парадокса Е.З. Демиденко [31], поскольку согласно ему результат идентификации не опирается на непустоту информационного множества задачи.
3.4. Аппроксимация информационного множества
Если информационное множество задачи идентификации непусто, то возникает задача наиболее адекватного его представления. Как уже упоминалось, информационное множество рассматриваемой задачи идентификации, т. е. множество решений ИСЛАУ, в случае своей непустоты является объединением не более чем 2П полиэдров. Длина его полного прямого описания, при котором скрупулёзно выписываются в каждом ортанте все уравнения ограничивающих Б (А, Ь) гиперплоскостей, может расти экспоненциально с размерностью интервальной системы п. По этой причине в случае, когда нас в качестве решения задачи идентификации не удовлетворяет единственная точка, разумно пойти по пути построения оценок для информационного множества. Тогда рассматривается задача приближённого описания (оценивания) множества Б (А, Ь) в соответствии со смыслом решаемой практической постановки.
В зависимости от практической задачи возможны различные способы оценивания множеств решений. К примеру, если в рассмотренной выше задаче идентификации объекта в условиях интервальной неопределённости мы собираемся использовать результаты для дальнейшей оптимизации результата проектирования какой-либо большой системы, то логично предоставить лицу, принимающему решение, по возможности, более широкий набор альтернатив, из которого может быть выбрано решение всей большой задачи. Тогда имеет смысл оценивать Б (А, Ь) с помощью подмножеств, т. е. таким образом, чтобы наверняка не захватить "лишние" векторы параметров, не имеющие отношения к идентифицируемому объекту. Описание различных методов внутреннего оценивания множеств решений интервальных систем линейных алгебраических уравнений дано в [3, 6, 20, 38].
С другой стороны, результаты идентификации могут использоваться для гарантированной оценки по полученной модели наибольших возможных отклонений выходов в процессе функционирования объекта. Тогда множество решений Б (А, Ь) следует оценивать "извне", с помощью объемлющих множеств, которые содержат все векторы резуль-
татов идентификации. Методы внешнего оценивания множеств решений интервальных систем линейных алгебраических уравнений излагаются, в частности, в работах [1-4, 8, 10, 37-39].
3.5. Выявление выбросов в данных
По определению выбросом в данных называется результат измерения, выделяющийся из общей выборки. Иными словами, выбросы в измерениях и наблюдениях — это такие результаты, которые заметно отклоняются от остальных членов выборки, в которой они присутствуют. Выбросы, как правило, являются следствиями аномальных ошибок измерений или наблюдений, результатом случайного и неконтролируемого воздействия на условия опыта и т. п., и потому их выявление является важной частью общего процесса обработки данных. Особую сложность процедуре определения выбросов придаёт тот факт, что выбросы могут быть неединственными, а их общее количество может быть неизвестным.
Характеристический признак выбросов состоит в том, что их удаление из выборки делает её существенно более согласованной, по крайней мере настолько, насколько это возможно при наличии других неустранимых ошибок. Отсюда следует универсальный метод выявления выбросов:
— из выборки удаляются заметно отличающиеся от остальных результаты измерений, исследуется согласование оставшейся части выборки;
— если данные в оставшейся части выборки имеют существенно лучшее согласование, чем в исходной, то удалённая часть выборки может считаться выбросом.
Однако трудоёмкость этой процедуры может быть очень большой. Так, чтобы выявить выброс из одного наблюдения в выборке из т наблюдений, нужно исследовать совместность т подвыборок размера (т — 1). Чтобы выявить два наблюдения-выброса в выборке из т наблюдений, требуется перебрать С! = т(т —1)/2 подвыборок из (т — 2) остающихся наблюдений. И так далее. В целом в этом процессе придётся исследовать
С! + ст + ст + ... подвыборок, где ст = -¡-Г/-гтг — число сочетаний из т по к.
к! (т — к)!
Размер подвыборок, которые имеет смысл удалять из исходной выборки в вышеописанной процедуре, логично брать не большим т/2, так как только тогда оставшаяся часть выборки образует "большинство наблюдений". Представление об общей трудоёмкости процесса даёт известное тождество С! + С! + С! + ... + С! = 2т.
В предлагаемом методе максимизации согласования информация о характере максимума распознающего функционала позволяет выявлять в выборке наблюдения, подозрительные на выбросы. Покажем, как это можно делать.
Пусть распознающий функционал
Шв (х) = шт < rad 5^ + (га^ а?) Х | —
¿=1
mid — ^^ (mid а?). ¿=1
достигает своего безусловного максимума, равного М, в точке т. Рассмотрим в этой точке значения образующих распознающего функционала, т. е. выражений
^¿(х) = rad + ^^ (rad а?) |х? | — ¿=1
mid 5^ — ^^ (mid а?) ¿=1
1, 2,
т.
Обозначим посредством Im множество индексов тех образующих которые принимают значение M при x = т; Im — это подмножество множества первых m натуральных чисел, которое, очевидно, непусто и содержит хотя бы один элемент — тот индекс, на котором достигается M. Тогда у остальных образующих, с номерами из дополнения { 1, 2,... , m} \ Im, значения в точке т больше, чем M. Поэтому
max Uss (x) = Uss (т)= min ^(т) < min ^(т) < max min ^(x).
x€Rn ie{1,2,...,m} ie{1,2,...,m}\JM x€Rn ie{1,2,...,m}\/M
Таким образом, если удалить из выборки все наблюдения, номера которых принадлежат множеству Im , то значение максимума распознающего функционала может только увеличиться. Это означает увеличение меры согласования оставшихся в выборке данных.
Как следствие, если в точке максимума распознающего функционала Uss значения некоторых, относительно немногих от общего числа, образующих функционала принимают существенно меньшие значения, чем остальные образующие, то это нужно расценивать как признак того, что соответствующие этим образующим наблюдения являются выбросами.
3.6. Практическая реализация
Практичность метода максимума согласования для решения задачи идентификации зависит главным образом от эффективности нахождения максимума распознающего функционала Uss. В самом общем случае, когда интервальная неопределённость присутствует на входе и выходе системы, распознающий функционал информационного множества задачи является негладкой и многоэкстремальной целевой функцией, так что её оптимизация весьма сложна. Для решения этой задачи применимы подходящие методы глобальной оптимизации, причём можно выбрать их адаптивными, т. е. подстраивающими своё выполнение под задачу, в отличие от пассивных подходов к исследованию разрешимости ИСЛАУ, упомянутых в разделе 2.1.
Важный частный случай задачи идентификации соответствует точному заданию входных величин а^, когда интервальная неопределённость присутствует лишь на выходах bj (рис. 7). В настоящее время этот случай наиболее детально разработан в теории нестатистического оценивания параметров (см., в частности, [21, 23, 24, 33, 40, 41]
Рис. 7. Для точных входных данных задача оценивания параметров сводится к проведению гиперплоскости через интервалы выходной величины
и цитированную в этих работах литературу). Отметим также, что при дополнительных статистических предположениях данный случай соответствует условиям применения традиционного регрессионного анализа, при которых получены наиболее сильные результаты об оптимальности популярного метода наименьших квадратов.
В случае, когда значения входных переменных точны, получаем интервальную линейную систему Ах = Ь с точечной матрицей А = (а.), в которой интервальность сосредоточена лишь в правой части. Как следствие, распознающий функционал принимает более простой вид
ивв (х, А, Ь)
Ш1п -ч rad Ь —
1< г<т
mid Ьг — ^^ <
.7=1
Теперь он ограничен сверху, так как все rad а. = 0 и, кроме того, в силу Предложения 2 является глобально вогнутым на всём пространстве Кга. Вместо многоэкстремальной конфигурации, изображённой на рис. 3, получаем унимодальный функционал, график которого выглядит примерно как на рис. 8. На этом рисунке показан график распознающего функционала множества решений для ИСЛАУ
( 3 -1\
-1 2
V
2 0
х1 Х2
Л—2, 2Л [0,1] [—1,0] V [0, 2] )
Хорошие свойства функционала Шя в случае точечной матрицы А позволяют использовать для отыскания его максимума развитые методы негладкой выпуклой оптимизации (см., например, [42, 43]). Для систем компьютерной математики йеНаЪ и МАТ1АВ свободно распространяемая программа lintregr, которая реализует эту версию метода максимума согласования, находится на сервере ИВТ СО РАН по адресу http://www.nsc.ru/interval/Programing. Она использует в качестве основы код
Рис. 8. График распознающего функционала для ИСЛАУ с точечной матрицей
ralgb5, написанный д.ф.-м.н. П.И. Стецюком (Институт кибернетики НАН Украины, Киев).
Другой возможный способ максимизации распознающего функционала в случае точечной матрицы — решить задачу линейного программирования по отысканию максимума u для пар (x,u) G Rn+1, х G Rn, u G R, принадлежащих выпуклому полиэдральному подграфику распознающего функционала Uss. Формулировка этой задачи имеет следующий вид:
!u < a(i)x — b., ,
(Л — i =1, 2,..., m, u < —awx + b.,
где a(i) означает i-ю строку матрицы A рассматриваемой ИСЛАУ. В целом для точных входных данных получаем эффективную методику обработки данных с интервальными неопределённостями в выходных переменных, которая является хорошей альтернативой методу наименьших квадратов. Она находит так называемую "чебышёвскую точку" системы линейных неравенств [44]. Данная методика очевидным образом обобщается на нелинейный случай (см. раздел 2.5), но требует привлечения существенно более сложных оптимизационных методов.
Если в общем случае l — общее количество входных переменных с неопределённостью, то, как показывает Предложение 2, распознающий функционал Uss имеет 21 областей вогнутости. Поэтому, когда l невелико, можно организовать нахождение глобального максимума распознающего функционала через перебор всех его областей вогнутости и максимизацию на каждой из них с помощью методов негладкой оптимизации, описанных в [42, 43]. Перспективной является также техника кодифференциального спуска [45, 46].
3.7. Сравнение с методом наименьших квадратов
В заключение этой части работы приведём сравнение метода максимума согласования с методом наименьших квадратов (МНК) на примере, особенностью которого является то, что оценка по МНК не лежит в информационном множестве задачи идентификации.
Предположим, что переменная y G R линейно зависит от переменной х G R по правилу
y = ах + ß.
Неизвестные коэффициенты а и ß надо определить по результатам измерений, данные которых приведены в следующей таблице:
Номер опыта 1 2 3
x 0 1 2
У 1 2 —0.5
Предполагается также, что в каждом опыте переменная х измеряется точно, тогда как для переменной у измерения дают только интервал, середина (центр) которого указана в таблице. Радиус интервала неопределённости во всех опытах равен единице, а настоящее значение у может быть любым числом из этого интервала. Требуется оценить информационное множество, т.е. множество пар значений а и в, согласующихся с данными измерений.
Введём обозначения: х — значение переменной х в г-м опыте, у^ — центр интервала для переменной у в г-м опыте.
Информационное множество в рассматриваемом примере описывается системой включений
На рис. 9 полосы пронумерованы, а информационное множество выделено заливкой. Это треугольник с вершинами (—1, 2), (-0.5,1.5) и (-0.75, 2).
Воспользуемся для оценивания коэффициентов а и в методом наименьших квадратов, предполагая, что неопределённость в измерениях выходной переменной у имеет теоретико-вероятностный характер [47]. Представленные в таблице выше значения у можно считать средними значениями (математическими ожиданиями), а интервалы неопределённости пусть содержат все или почти все возможные значения этой величины. Например, если ошибки в у подчинены нормальному закону распределения с дисперсией а, то интервалы ошибок [-1,1] для у могут быть хорошо известными инженерам интервалами "плюс-минус За" или даже "плюс-минус 6а", которые обеспечивают достоверность на уровне 99.73 и 99.99966 % соответственно.
(24)
и представляет собой пересечение трёх полос:
(I) в е [0, 2], (II) в е-а + [1, 3], (III) в е-2а + [-1.5, 0.5].
а
Рис. 9. В пространстве переменных а и в оценка по МНК (звёздочка) не лежит в информационном множестве (затенённый треугольник)
Обозначим через а* и в* оценки по МНК для а и в. Величины а* и в* определяются, как известно, из нормальной системы уравнений
112) (0 ■) о=с 1 э (—0.) ■ (25)
которая получается из средней системы уравнений для (24) домножением слева на транспонированную матрицу коэффициентов. Как следствие (25), имеем
5 3) /аЛ _ ( 1 3 у \ву = \2.5
1 /5 3\ /5 3)-1 1 ( 3 —3
deЦ з I = 6, обратная матрица — (33) = 6 I -3 5
так что искомая оценка
/а*) =1 ( 3 — 3) ( 1 ) = 1 /—■4.5) = ( —3/4) = /—0.75 ) = 6 ^—3 5 ) \2.5y = 6 V 9.5 ) = \,19/12у V 1.58(3))
Точка с координатами (а*, в*) показана на рис. 9 звёздочкой и лежит на заметном удалении от информационного множества задачи. Сравнение оценки по МНК с информационным множеством в координатах (х,у) представлено на рис. 10.
Результаты оценивания по предложенному в работе интервальному методу максимизации согласования таковы:
шахивв = 0.125, аггшахивв = ( ^^
1.875
Рис. 10. В пространстве переменных х и у прямая у = а*х + в*, определяемая с помощью МНК (пунктир), не лежит в множестве всех прямых, проходящих через интервалы измерений
Как видим из рис. 9, полученная оценка (жирная точка) находится в информационном множестве. На рис. 10 этой точке соответствует выделенная штрихпунктиром регрессионная прямая, расположенная примерно посредине всех возможных прямых, согласованных с данными.
Отметим, что использование взвешенного метода наименьших квадратов (см., например, [48]) с подходящими весами различных наблюдений позволяет получить результат, который лежит в информационном множестве задачи. Но выбор этих нужных весов является самостоятельной нетривиальной задачей.
Выводы
Введением распознающего функционала множества решений задача исследования разрешимости интервальной системы линейных алгебраических уравнений сводится к удобной аналитической форме, позволяющей более детально исследовать и корректировать исходную систему. В работе предложен распознающий функционал Uss, который имеет определённые преимущества перед функционалом Uni, предложенным ранее.
Метод максимума согласования — перспективный метод идентификации параметров систем и обработки данных с интервальными неопределённостями, основанный на максимизации распознающего функционала множества решений (информационного множества задачи). Он может служить хорошей альтернативой традиционным методам регрессионного анализа, использующим теоретико-вероятностные модели ошибок данных.
Список литературы
[1] Алефельд Г., Херцбергер Ю. Введение в интервальные вычисления. М.: Мир, 1987.
[2] Калмыков С.А., Шокин Ю.И., Юлдашев З.Х. Методы интервального анализа. Новосибирск: Наука, 1986.
[3] Шарый С.П. Конечномерный интервальный анализ. Новосибирск, 2013. Электронная книга, доступная на http://www.nsc.ru/interval/Library/InteBooks
[4] Neumaier A. Interval Methods for Systems of Equations. Cambridge: Cambridge Univ. Press, 1990.
[5] Shary S.P. A new approach to the analysis of static systems under interval uncertainty // Scientific Computing and Validated Numerics. Proc. of the Intern. Symp. on Sci. Comput., Computer Arithmetic, and Validated Numerics SCAN'95 held in Wuppertal. Germany, September 26-29, 1995 / Eds. G. Alefeld, A. Frommer and B. Lang. Berlin: Akademie Verlag, 1996. P. 118-132.
[6] Шарый С.П. Алгебраический подход к анализу линейных статических систем с интервальной неопределённостью // Изв. РАН. Теория и системы управления. 1997. №3. С. 51-61.
[7] Shary S.P. A new technique in systems analysis under interval uncertainty and ambiguity // Reliable Comput. 2002. Vol. 8, No. 5. P. 321-418. (Электронная версия статьи доступна на http://www.nsc.ru/interval/shary/Papers/ANewTech.pdf)
[8] Фидлер М., Недома Й., Рамик Я. и др. Задачи линейной оптимизации с неточными данными. Москва; Ижевск: Изд-во РХД, 2008.
[9] Kearfott R.B., Nakao M., Neumaier A. et al. Standardized notation in interval analysis // Вычисл. технологии. 2010. Т. 15, № 1. С. 7-13.
[10] Moore R.E., Kearfott R.B., Cloud M.J. Introduction to Interval Analysis. Philadelphia: SIAM, 2009.
[11] Ерёмин И.И. Противоречивые модели оптимального планирования. М.: Наука, 1988.
[12] Лакеев A.B., Носков С.И. О множестве решений линейного уравнения с интервально заданными оператором и правой частью // Сибирский матем. журн. 1994. Т. 35, № 5. С. 1074-1084.
[13] Kreinovich V., Lakeyev A., Röhn J., Kahl P. Computational Complexity and Feasibility of Data Processing and Interval Computations. Dordrecht: Kluwer, 1998.
[14] Шарый С.П. О характеризации объединённого множества решений интервальной линейной алгебраической системы. Красноярск, 1990. 20 с. Деп. в ВИНИТИ, № 72б-В91.
[15] Shary S.P. Solving the linear interval tolerance problem // Math. and Comput. in Simulation. 1995. Vol. 39. P. 53-85.
[16] Шарый С.П. Решение интервальной линейной задачи о допусках // Автоматика и телемеханика. 2004. № 10. С. 147-1б2.
[17] Шарый С.П. Разрешимость интервальных линейных уравнений и анализ данных с неопределённостями // Там же. 2012. № 2. С. 111-125.
[18] Beeck H. Über die Struktur und Abschatzungen der Lösungsmenge von linearen Gleichungssystemen mit Intervallkoeffizienten // Computing. 1972. Vol. 10. P. 231-244.
[19] Рокафеллар Р. Выпуклый анализ. М.: Мир, 1973.
[20] Шарый С.П. Ещё раз о внутреннем оценивании множеств решений интервальных линейных систем // Вычисл. технологии. 2003. Т. 8, Спец. выпуск "Труды Совещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям". С. 14б-1б0.
[21] Вощинин А.П., Бочков А.Ф., Сотиров Г.Р. Метод анализа данных при интервальной нестатистической ошибке // Заводская лаборатория. 1990. Т. 5б, № 7. С. 7б-81.
[22] Вощинин А.П. Интервальный анализ данных: Развитие и перспективы // Там же. 2002. Т. б8, № 1. С. 118-12б.
[23] Вощинин А.П. Задачи анализа с неопределёнными данными — интервальность и/или случайность? // Междунар. конф. по вычислительной математике МКВМ-2004. Рабочие совещания / Под ред. Ю.И. Шокина, А.М. Федотова, С.П. Ковалёва, Ю.И. Молоро-дова, А.Л. Семёнова, С.П. Шарого. Новосибирск: ИВМиМГ СО РАН, 2004. С. 147-158. (Электронный адрес статьи http://www.nsc.ru/interval/Conferences/IMRO_G4/Voschi-nin.pdf)
[24] Жилин С.И. Нестатистические модели и методы построения и анализа зависимостей. Дис. ... канд. физ.-мат. наук. Барнаул: АлтГУ, 2004. 119 с. (Электронная версия доступна по адресу http://www.nsc.ru/interval/Library/ApplDiss/Zhilin.pdf)
[25] Носков С.И. Технология моделирования объектов с нестабильным функционированием и неопределенностью в данных. Иркутск: Облинформпечать, 199б.
[26] Подружко A.A., Подружко А.С. Интервальное представление полиномиальных регрессий. М.: Эдиториал УРСС, 2003.
[27] Nguyen H.T., Kreinovich V., Wu B., Xiang G. Computing Statistics Under Interval and Fuzzy Uncertainty. Berlin; Heidelberg: Springer-Verlag, 2012.
[28] Канторович Л.В. О некоторых новых подходах к вычислительным методам и обработке наблюдений // Сибирский матем. журн. 19б2. Т. 3, № 5. С. 701-709.
[29] Куржанский A^. Задача идентификации — теория гарантированных оценок // Автоматика и телемеханика. 1991. № 4. С. 3-2б.
[30] Кумков С.И. Обработка экспериментальных данных ионной проводимости расплавленного электролита методами интервального анализа // Расплавы. 2010. № 3. С. 79-89.
[31] Демиденко Е.З. Комментарий II к статье А.П. Вощинина, А.Ф. Бочкова и Г.Р. Соти-рова "Метод анализа данных при интервальной нестатистической ошибке" // Заводская лаборатория. 1990. Т. 56, № 7. С. 83-84.
[32] Хлебников А.И. О проблемах использования метода центра неопределённости для обработки экспериментальных данных // Вычисл. технологии. 1999. Т. 4, № 4. С. 80-81.
[33] Оскорбин Н.М., Максимов А.В., Жилин С.И. Построение и анализ эмпирических зависимостей методом центра неопределённостей // Изв. Алтайского гос. ун-та. 1998. № 1. С. 35-38.
[34] Zhilin S.I. On fitting empirical data under interval error // Reliable Comput. 2005. Vol. 11, No. 5. P. 433-442.
[35] Zhilin S.I. Simple method for outlier detection in fitting experimental data under interval error // Chemometrics and Intellectual Laboratory Systems. 2007. Vol. 88, No. 1. P. 60-68.
[36] Суханов В.А. Исследование эмпирических зависимостей: Нестатистический подход. Барнаул: Изд-во Алтайского гос. ун-та, 2007.
[37] Поляк Б.Т., Назин С.А. Оценивание параметров в линейных многомерных системах с интервальной неопределённостью // Проблемы управления и информатики. 2006. № 1-2. С. 103-115.
[38] Жолен Л., Кифер М., Дидри О., Вальтер Э. Прикладной интервальный анализ. Москва; Ижевск: Изд-во РХД, 2007.
[39] Gay D.M. Solving interval linear equations // SIAM J. Numer. Anal. 1982. Vol. 19, No. 4. P. 858-870.
[40] Померанцев А.Л., Родионова О.Е. Построение многомерной градуировки методом простого интервального оценивания // Журн. аналит. химии. 2006. Т. 61. С. 1032-1047.
[41] СпивАк С.И., ИсмАгиловА А.С. Математические модели химической кинетики. Нефтекамск: РИЦ БашГУ, 2010.
[42] Шор Н.З., Журбенко Н.Г. Метод минимизации, использующий операцию растяжения пространства в направлении разности двух последовательных градиентов // Кибернетика. 1971. № 3. С. 51-59.
[43] Шор Н.З., Стеценко С.И. Квадратичные экстремальные задачи и недифференцируе-мая оптимизация. Киев: Наук. думка, 1989.
[44] Зуховицкий С.И., Авдеева Л.И. Линейное и выпуклое программирование. М.: Наука, 1967.
[45] Angelov T. Solvability of systems of interval linear equations via the codifferential descent method // 15th GAMM-IMACS Intern. Symp. on Sci. Comput., Comput. Arithmetics and Verified Numerics SCAN'2012. Novosibirsk, ICT SB RAS, Russia, September 23-29, 2012. Book of Abstr. P. 13-14.
[46] Ангелов Т.А. О вычислении кодифференциалов // Вычисл. методы и программирование. 2013. Т. 14. С. 113-122.
[47] Линник Ю.В. Метод наименьших квадратов и основы теории обработки наблюдений. 2-е изд. М.: ГИФМЛ, 1962.
[48] Зоркальцев В.И. Метод наименьших квадратов: Геометрические свойства, альтернативные подходы, приложения. Новосибирск: Наука, 1995.
Поступила в 'редакцию 21 января 2013 г., с доработки — 27 марта 2013 г.
