CC BY
30Метод нечеткой линеаризации для численного решения алгебраических и дифференциальных уравнений
Евдокимов А.В. feval@applmech.mipt.ru )
Московский Физико-Технический Институт (Государственный университет)
1. Постановка задачи
Нечеткие методы применяются при математических расчетах в слабо формализуемых предметных областях, когда нельзя пользоваться стохастическими методами (собранными, например, в [9]) ввиду неприменимости понятия случайной функции и отсутствия информации о распределениях вероятностей параметров. Нечеткие числа могут быть заданы в различных формах (иметь различные функции принадлежности). В частности, к ним относятся АК-числа [3], гауссовские числа (нормально распределенные случайные величины с заданной дисперсией), числа, представленные в виде набора интервалов [2] и т. д.
Для каждой из этих форм существует несколько вариантов алгебры, которые напрямую используются в численных методах решения алгебраических уравнений [2,3,10] и задач оптимизации [7]; при этом часто даже не требуется модификация алгоритма, предназначенного для решения данной задачи в вещественных числах. Однако существующие алгебры либо определяют функцию принадлежности результата операции на максимально широком носителе (например, стандартный интервальный анализ), что обеспечивает математическую строгость (но при решении практических задач дает завышенную степень неопределенности) [1], либо необоснованно уменьшают неопределенность результата искусственными приемами. К алгебраическим подходам второго типа, которые более применимы на практике, но могут давать некорректные результаты, относятся многочисленные модификации интервальной алгебры («нестандартные вычитание и деление», мнимые интервалы, квазилинейное пространство и т. п.) [1,2,4,6], а также приближенная арифметика АК-чисел [3]. Найти удовлетворительный компромисс между строгостью, практичностью и универсальностью нечетких расчетов пока удается только в простых частных случаях — в частности, в безытерационных методах решения линейных алгебраических уравнений с некоторыми ограничениями (например, интервальный метод прогонки) [5,6]. Несмотря на то, что численные методы для дифференциальных уравнений выглядят как алгебраические соотношения, задачу использования при их расчете какой-либо алгебры нечетких чисел (без модификации четкого численного метода) пока нельзя назвать решенной [8].
Поэтому сейчас более перспективным считается альтернативное направление развития нечеткой математики (уже сугубо численное), в котором функции принадлежности нечетких результатов (как правило, разложенные на интервалы) восстанавливаются на основе многократных решений четких задач с комбинированием значений параметров [11,13]. Этот подход идейно близок численном методам типа Монте-Карло [9] для решения стохастических уравнений, обладает по сравнению с ними более широкой областью применения, но требует еще более значительных вычислительных ресурсов. Он предполагает одновременное хранение и обработку очень большого числа «копий» переменных (порядка М , где К - число нечетких параметров); попытки уменьшить это число (за счет М, с последующей интерполяцией) математически корректны лишь при малых К. Несомненным преимуществом таких методов является лишь выгодность применения параллельных вычислений.
В обоих указанных подходах либо используется фиксированное представление нечетких чисел, либо произвольное представление с той или иной степенью точности аппроксимируется набором интервалов. В то же время, при внедрении нечеткого анализа в
программные пакеты моделирования крайне желательным свойством является независимость расчетов от конкретного способа формализации нечеткости (выбор которого зависит от задачи, неизвестной при разработке пакета). Поэтому актуальной является задача такой формулировки методов решения дифференциальных и нелинейных алгебраических уравнений, которая бы не зависела от формы нечетких чисел.
Предлагаемый ниже метод является инвариантным относительно формы чисел, позволяет по сравнению с существующими чисто алгебраическими методами более точно решать относительно сложные нечеткие задачи — этим он похож на методы сведения к системе четких задач, — но по сравнению с ними обладает большей универсальностью и производительностью (свойствами, присущими алгебраическим методам): см. таблицу 1. Основная преодолеваемая в данном методе проблема алгебраических методов состоит в том, что арифметические операции в численных методах всегда проводятся не только над независимыми нечеткими числами, но и над числами, являющимися результатами операций над одними и теми же нечеткими исходными данными. Стандартная нечеткая алгебра дает заниженную характерную погрешность результирующего нечеткого числа в случае сложения и умножения зависимых чисел, и завышенную — в случае вычитания и деления. Нестандартные же алгебры (например, алгебра мнимых интервалов [1]), напротив, независимые нечеткие числа обрабатывают как зависимые.
Таблица 1. Сравнение предлагаемого метода нечеткой линеаризации с существующими подходами к численному решению нечетких уравнений
Подход Сравнение метода нечеткой линеаризации с подходом
Общее Преимущества Недостатки
I. Сведение Сопоставимая 1. Отсутствие а) привязки к 1. Слабая
к системе степень классу задач, б) ограничений теоретическая
четких задач неопределенности на число нечетких база.
результатов параметров, в) необходимости аналитических выкладок и интервального представления. 2. Высокая производительность 2. Более приближенное описание нелинейностей
II. На основе 1. Высокая 1. Расчет не максимально Отсутствие
строгой алгебры универсальность. возможной, а «характерной» гарантированной
2. Использование неопределенности. принадлежности
стандартных 2. Возможность деления на к множеству-
(«четких») число, интервальное результату
численных представление которого
методов содержит ноль
III. На основе —//—. 1. Автоматическое Большие
«исскуственной» Уменьшение различение зависимых и затраты памяти
алгебры неопределенности независимых нечетких и машинного
учетом операций с чисел. 2. Использование времени
зависимыми любого представления
числами чисел. 3. см. 11-2
Таким образом, целью данной работы является построение приближенного метода проведения алгебраических операций над нечеткими числами с учетом зависимости чисел-операндов, применимого (в сочетании с произвольными четкими численными методами) для решения нечетких алгебраических и дифференциальных задач. Метод ориентирован на преодоление проблемы экспоненциального роста вычислительных затрат с числом нечетких
*
Некоторые преимущества относятся только к прямому методу линеаризации.
параметров, которая характерна для современных неалгебраических методов. Метод имеет сходную формулировку и дает близкие результаты для различных форм нечетких чисел. На примере гауссовских и интервальных чисел эффективность метода в работе показывается в применении к нечеткой алгебраической модели сердца и к нескольким дифференциальным задачам из различных областей.
2. Метод нечеткой линеаризации для расчета нечетких уравнений
2.1. Метод линеаризации как метод учета зависимостей чисел
Идея преодоления проблемы операций над зависимыми нечеткими числами состоит в том, что число должно хранить не только свое текущее значение (включающее некоторым образом формализованную погрешность), но и информацию о том, из каких исходных данных и как это число было получено. Это дополнительно дает очень полезную в приложениях возможность анализа результатов численных экспериментов на предмет того, каким образом сказались на них заданные исходные данные. Однако буквальная реализация указанной идеи приводит к неадекватным затратам памяти и времени расчетов: фактически, нечеткое число превращается из значения в формулу зависимости значения от исходных данных; причем с каждым параметром этой формулы при любой арифметической операции должны производиться сложные вычисления.
В целях экономии ресурсов любое нечеткое число (независимо от способа его формализации) предлагается представлять в виде линейной комбинации по нечетким числам - исходным данным, что выражается формулой
к,
X = х0, + Е с,у %у - х0, + Е сцк % ]к , (1)
МЛ{ к=1
где х0, — часть числа х,, не зависящая от исходных нечетких чисел Ъу. В числе (х,) должны храниться лишь ссылки на эти числа (или их идентификаторы у, в зависимости от языка реализации), и скалярные коэффициенты Су при них. Ниже такая конструкция называется «линеаризованной историей» числа. В отличие от отмеченного выше полного варианта хранения информации об истории, текущее значение нечеткого числа (формализованное одним из указанных во введении способов) при использовании линеаризованной истории не теряет смысла и также должно храниться в числе.
Алгоритм арифметической операции над числами х1 и х2 (,=1,2) таков:
1. Определяется набор исходных данных {}, которые должны входить в линеаризованную историю результата операции. В простейшем случае это делается объединением множеств {/}1 и {/}2, хотя число элементов множества {} может быть уменьшено (см. 2.3).
2. Для каждого 7'е{/}1П{/}2 рассчитываются коэффициенты линейной комбинации (для остальных у из {/}1 или {у}2 расчет тривиален). В случае сложения/вычитания это можно сделать точно:
х = х1 ± х2 ^ с± = с1 у ± с2у, (2а)
в случае умножения/деления — лишь приближенно (заменяя по очереди каждый из нечетких операндов х, на его среднее скалярное значение а,):
х = х1х2 ^с* « ЦС1 уа2 + (1 - ^)с2уа1, (2б)
где вес q является одинаковым для всех у, в простейшем случае равен 1/2 и может зависеть от всех произведений суа3-,. Выбор весовой функции q влияет на погрешность метода. В случае деления можно использовать дополнительное (по сравнению с умножением) приближение 1/х~х/а2 (которое часто встречается в нечеткой алгебре):
— * / 2
х = х1/х2 ^ су « су а2 , (2в)
однако вместо этого рекомендуется рассматривать 1/х как элементарную функцию (см. формулу (6)), вследствие чего
— * / 2
х = х1/х2 ^ с/ « - с / а,2 . (2в')
3. Вычисляется величина поправки к погрешности, обусловленной наличием в линеаризованных историях х1 и х2 одних и тех же чисел Ъ/. Данный шаг алгоритма является единственным, который требует проведения различных аналитических выкладок для разных форм нечетких чисел. В случае гауссовских чисел эти выкладки (см. ниже) дают для сложения/вычитания
Да2 = ±2^5+ , 5+ = С1 с/а2 (3а)
}
для умножения —
Да2 = , 5* = (с1 }С2/ У (а/а} У (3б)
}
а для деления —
Да2 =-41, (3в)
а2 /
где о2/ — дисперсии исходных нечетких чисел, Да2 — (аддитивная) поправка к дисперсии о2, вычисляемой на шаге 4. 4. По правилам соответствующей нечеткой алгебры проводится обычная арифметическая операция (с числами х1 и х2 как с независимыми). Например, для гауссовского числа на этом шаге вычисляется как среднее значение а, так и его дисперсия о2; для сложения/вычитания независимых чисел х1 и х2 имеем
о2 = о12 + о22, (4а)
для умножения/деления —
о2 = а2(о\/а\ + 022/а22). (4б)
Итоговая погрешность (в данном случае — дисперсия) рассчитывается с учетом найденной на шаге 3 поправки.
Правила арифметических операций между нечетким числом и скаляром достаточно очевидны: при сложении/вычитании коэффициенты линеаризованной истории нечеткого числа не изменяются, а при умножении/делении они умножаются/делятся на этот скаляр.
Для доказательства формул коррекции погрешности (3) обозначим он2 дисперсию, получаемую в предположении о независимости операндов, а £н — слагаемые, соответствующие исходным числам входящим только в один из операндов: дае{/}, {/}1П{/}2. Тогда в случае сложения/вычитания
а2 =а2н + Да2 =^(С/а / )2 + Sн, )
аН = /а/) +^(с2/а/) + Sн, следовательно, так как С/ = с1/- ± с2/, 1 )
Да2 = ^ [(с1 / ± с2/ ) - с2/ - с|/ ]а2 = ±2^с1 /с2/а2 , что и требовалось доказать.
1 )
Доказательство в случае умножения базируется на формуле дисперсии произведения 2 2 2 2 2 2 2 независимых (о = а (о1 /а1 + о2 /а2 )) и на формуле дисперсии функции (о(х ) = |2х|о(х)):
а2 = Д2 )У + ^ = 4£(/ £/а /)2 + ^,
= Еа2 (2 &)+бн = е]] ] )2 НД2-+а2Д2)+ ъ ,
] ]
Аа2 = а2 -а2 = ^с^ у а ])2 - Т^] % - ])2 •
] ] ]
Формула для деления (3в) не следует из формулы для умножения (хотя они отличаются только множителем -а24). Ее вывод здесь не приводится по причине своей громоздкости.
2.2. Прямой вариант метода нечеткой линеаризации
Следует заметить, что в стандартной нечеткой алгебре (гауссовских и большинства других нечетких чисел) отсутствует свойство дистрибутивности погрешности (ширины интервала и т. п.). В частности, если х = а-(Ъ+с), а у = а-Ъ+а-с, то алгебра независимых гауссовских чисел (см. формулы 4) дает (ох2-оу2) = 2-Ъ-с-оа2. Поправка из формулы (3а), учитывающая, что оба числа а в формуле для у — это одно и то же число, возвращает свойство дистрибутивности погрешности. Однако это верно только в том случае, если нечеткие числа, содержащие в своей линеаризованной истории одно и то же исходное нечеткое число, входят в формулы аддитивно (не перемножаясь). В противном случае даже предлагаемая поправка к погрешности (см. формулу 3б) не обеспечивает свойство дистрибутивности. Точнее, это свойство можно обеспечить за счет выбора весового коэффициента q, но далеко не всегда и лишь ценой весьма ресурсоемких расчетов. Нарушение дистрибутивности допустимо в безытерационных алгоритмах решения алгебраических уравнений, однако в более сложных задачах оно приводит к качественно неверному решению, — например, к экспоненциальному росту погрешности (со временем или с номером итерации). Если в формулах итерационных методов встречается перемножение зависимых нечетких чисел, то вместо шагов 3-4 описанного выше алгоритма предлагается явно применять правила соответствующей алгебры по отношению к исходным нечетким числам. В случае гауссовских чисел это дает:
а2 (х1 + Х2 )=Е с 2 а2, (5а)
)
а2 (х1 х2 ) = ЕЕ ( с2 ]2 } а2 ( £ ]2 )+ х01а 2 + х02а2 , где (5б)
]1 ]2
2 / \ 2 2 9 9
а ((1 £] ) = а22 + а21 для
а2 ((1 £] )= 4а2а2 для] 1=j9=j•
Эти формулы получаются напрямую путем сложения (5а) и умножения (5б) чисел х1 и х2 в представлении (1). Хотя необходимость в таком прямом методе вызвана свойствами операции умножения, наряду с (5 б) необходимо также использовать формулу для сложения (5а) вместо (3а,4а). Сложение (5а) требует даже меньше машинных операций, чем (3а), однако умножение в прямом методе является гораздо менее экономичным: число машинных операций становится пропорциональным К2. По этой причине для грубых оценок целесообразно производить умножение нечетких чисел не по старым, а по новым коэффициентам линейной комбинации, заменяя
К2
умножений К сложениями:
а 2 (х1 Х2 )=Е с? а 2. (5б')
)
Преимуществом формул (5) по сравнению с (3-4) является отсутствие необходимости в аналитических выкладках (см. раздел 2.1), которые для некоторых типов нечетких чисел являются очень громоздкими, а для произвольных функций принадлежности без интервального представления вообще невозможны. Кроме того, в прямом варианте метода достаточно решить задачу один раз при каких-то одних значениях параметров, заданных в какой-либо одной форме (можно даже при четких значениях), после чего можно много раз
подставлять в найденные линейные комбинации не только другие значения нечеткости, но и нечеткие числа (с тем же средним), представленные в другой форме. Другими словами, каждый следующий расчет одной нечеткой задачи может проводиться без пересчета коэффициентов с] и даже без использования базовых вычислительных алгоритмов. Это существенно экономит машинное время в нечетких вычислительных экспериментах.
Предлагаемое правило изменения коэффициентов при вычислении произвольной элементарной функции над нечетким числом у = /(х) сводится к их умножению на скалярное значение производной этой функции
су] * I '(аУх] . (6)
Обосновать это приближенное правило можно исходя как из предположения о малости вклада каждого ' в х, так и из линейной экстраполяционной формулы (разложения в ряд Тейлора до первого члена) для 1(х) в целом. Для того чтобы сумма линейной комбинации вычисленного таким образом результата функции была равна скалярному значению результата, необходимо также положить
у0 = 1 (а) су/а] . (6)
)
Если производная /'(а) определяется численно, то формулы (6,6') можно рассматривать как интерполяцию значения функции нечеткого аргумента по близким (четким) соседним точкам на кривой 1(а). Данный способ вычисления коэффициентов комбинации, не зависящий от функции 1(х), позволяет определить и само нечеткое значение функции, используя только ее аналог для вещественных чисел 1(а). Для этого достаточно восстановить нечеткость значения функции по найденной линейной комбинации, пользуясь стандартной формулой сложения независимых нечетких чисел (см. формулу 5а для гауссовских чисел).
Таким образом, прямой вариант метода нечеткой линеаризации является простой надстройкой над любой алгеброй независимых нечетких чисел (его реализация не требует изменения кода алгебраических операций), а также предъявляет очень малые требования к этой алгебре (в частности, не требует специальных алгоритмов вычисления нечетких функций).
2.3. Экономичность метода линеаризации
Количество машинных операций, необходимых для одной операции над нечеткими числами, пропорционально числу элементов К во множестве {/}1П{/}2, т. е. длине линеаризованной истории результата операции: К<К1+К2. Коэффициент пропорциональности для случая гауссовских чисел меняется в пределах от 4 до 13 (в зависимости от операции). Это является преимуществом рассмотренного метода, однако следует учесть, что в случае мультипликативного вхождения зависимых нечетких чисел в итерационные формулы целесообразно использовать прямой метод расчета погрешности умножения (5б), количество операций в котором пропорционально К1К2 (имеет вид 5К+(6^8)К1К2). При этом сложение в прямом методе занимает всего 3К операций, а вычисление нечеткой функции методом линеаризации, — 3К операций плюс затраты на однократное вычисление соответствующей скалярной функции и ее производной.
Преимущество прямого метода с точки зрения экономичности проявляется тогда, когда необходимо провести серию однотипных расчетов с одним и тем же скалярным значением нечетких параметров. В этом случае достаточно один раз провести полный расчет с указанными выше затратами, а затем для каждого набора исходных нечетких данных (которые могут отличаться, в том числе, и способом формализации нечеткости) затратить по 1К четких умножений и 1К четких сложений на каждое результирующее значение (для гауссовских чисел; здесь К — длина линеаризованной истории). При этом в случае дифференциального уравнения число результирующих значений может быть существенно меньше числа шагов по времени, поскольку для наглядного представления решения достаточно использовать лишь часть моментов времени.
В задачах с большим количеством исходных нечетких данных длина линейной комбинации К достаточно велика (числа рано или поздно будут зависеть от всех исходных данных). Поэтому вторым приемом, использующимся вместе с линеаризацией для уменьшения ресурсоемкости метода в итерационных задачах, является оставление в линейной комбинации только тех членов, которые наиболее существенны с точки зрения расчета погрешности. Существенность вклада исходного числа в погрешность определяется некоторой величиной Sj, которая для всех операций (см. формулы 3 и 5) пропорциональна модулю коэффициента с] при числе '■. Сокращение линеаризованной истории возможно за счет удаления из нее членов с наименьшими значениями этих величин (таких, что их относительный вклад в погрешность не больше заданного).
3. Результаты решения нечетких дифференциальных уравнений
Ниже с целью исследования свойств предложенного метода и с целью оценки его реальной производительности приводятся результаты расчетов с его помощью нескольких нечетких обыкновенных дифференциальных уравнений. Как элемент этих расчетов, с помощью метода также решались нечеткие алгебраические уравнения, однако сами по себе алгебраические уравнения будут рассмотрены в разделе 4.
3.1. Колебательные системы. Влияние четкого численного метода
Рассмотрим колебательные уравнения 2-го порядка: нелинейное уравнение Релея
ёх/ё = а ■ у + Ь ■ (х - х3 ) ёу/&
= - х
и его линейный аналог - простой осциллятор ёх/ё = а ■ у; ёу/ё = -х .
Соответствующие задачи Коши имеют три общих параметра: множитель а = 0.2±0.1 и начальные условия х0 = 2±0.5, у0 = 1±0.5. На рис. 1 показаны среднеквадратичные отклонения решений обоих уравнений, полученные при расчете на отрезке [0;80] с постоянным шагом 0.1 явным методом Рунге-Кутты 4-го порядка. На каждом графике приведены результаты пяти вариантов расчета, отличающихся набором нечетких параметров.
О 20 40 60 80 0 20 40 60 80
А. Линейный осциллятор Б. Уравнение Релея (Ь = аср)
Рис. 1. Погрешность (о) х-компонент решений колебательных дифференциальных уравнений. Цветом обозначены учитываемые нечеткие параметры: красный - х0, синий -у0, лиловый - (х0, у0), голубой - а, черный - (а, х0, у0).
На рис. 2 для линейного осциллятора результаты расчета погрешности сравниваются с модулем разности четких решений, полученных при значениях параметров, которые отличаются от указанных выше значений на ог- (в обе стороны). Тот факт, что кривая погрешности нечеткого решения во всех случаях (с 1-м, 2-мя и 3-мя нечеткими параметрами) ограничена кривыми разностей четких решений (либо совпадает с ними), показывает, что
предлагаемый подход дает результат, близкий к стандартной методике учета нечеткости за счет многократного расчета четких задач.
Из рис. 1-2 видно, что в случае нелинейного осциллятора наблюдается затухание колебаний погрешности с малой частотой, обусловленной изменениями разности фаз между решениями четких задач с разными собственными числами (что обусловлено «диссипативными» свойствами операции умножения с q = 1/2). Высокочастотные колебания расчетной погрешности даже в сильно нелинейных задачах (с Ъ > 1) совпадают с оценками погрешности по частоте и амплитуде, хотя форма колебаний может существенно отличаться.
Затраты машинного времени при использовании метода линеаризации оказались значительно меньше затрат на многократное решение четкой задачи (условия измерения и сравнения затрат см. в разделе 3.2). В частности, расчет уравнения Релея с одним нечетким а требует 0.71 с, с 2-мя нечеткими (х0у0) — 0.81 с, с 3-мя нечеткими (а,х0,у0) — 1.01 с; в то время как один четкий расчет занимает 0.25 с, а его 5, 25 и 125-кратные повторения занимают, соответственно, 0.86 с, 3.8 с и 19 с, что в 1.21, 4.7, и 19 раз больше.
о.6
0.45
"1 I Г
Л * * !?. Й *!Л Д *
0.15
О 20 40 60 80 0 20 40 60 30 " 0 20 40 60 80
А. Нечеткое значение а Б. Нечеткие значения (х0у0) В. Нечеткие значения (а,х0у0)
Рис. 2. Сравнение погрешности х-компоненты решения уравнения осциллятора (черная сплошная линия) с модулем разности решений, полученных при вариациях параметров
Л
0.01
0.005 -
-0.005
-0.01
0 20 40 60 А. Разность решений
_1_
_1_
_1_
0 20 40 60 80 Б. Разность погрешностей
20 40 60 80 В. Погрешности
Рис. 3. 7-компонента решения линейного осциллятора, полученного с помощью разных методов: 3-го порядка (красная кривая), 2-го (синяя) и 1-го (голубая). В качестве вычитаемого в разностях взято решение методом 4-го порядка (черная кривая на рис. В). На рис. А кривые 1-го и 2-го порядка для сопоставимости с 3-м порядком умножены на 0.02
0.02
0.01 -
-0.01 -
-0.02
0.02
0.01
0 20 40 60 80 А. Разность средних решений
-0.01 -
-0.02
0 20 40 60 80 Б. Разность погрешностей решений
Рис. 4. Разности у-компонент решений уравнения Релея (Ь = 0.1), полученных с помощью разных методов: 3-го порядка (красная кривая), 2-го (синяя) и 1-го (голубая). В качестве вычитаемого взято решение методом 4-го порядка.
Кривые 1-го и 2-го порядка для сопоставимости с 3-м порядком умножены на 0.02
Следует обратить внимание на то, что при использованном (относительно грубом) шаге по времени, равном 0.1, скалярное решение уравнения Релея (при Ь > 0.1) начинает расходиться к моменту времени около 50. Однако это не мешает корректно рассчитывать нечеткость решения, так как расходимость проявляется только в свободном члене линейной комбинации, но не в ее коэффициентах.
Выбор численного метода может существенно влиять на расчет погрешности решения даже в том случае, когда сочетания параметров задачи и шага интегрирования таково, что само решение (в четком смысле) слабо зависит от численного метода. Для иллюстрации этого факта уравнения линейного и нелинейного осциллятора решались (помимо рассмотренного выше явного метода типа Рунге-Кутты 4-го порядка) двумя неявными (точнее, полуявными) методами типа Рунге-Кутты: ^-устойчивым методом Розенброка 3-го порядка и монотонным ^-,^-устойчивым методом 2-го порядка. Для решения системы алгебраических уравнений на каждом временном шаге в линейном случае использовался прямой метод Гаусса, в нелинейном - метод Ньютона (с обращением матрицы Якоби тем же методом Гаусса). Сама матрица Якоби (матрица производных правой части системы уравнений) в разных вариантах вычислялась как по аналитической формуле, так и численно (простейшей аппроксимацией первого порядка точности), однако это практически не повлияло на результаты.
На рис. 3 и 4 показана разница между решениями колебательных уравнений явным методом 4-го порядка (см. выше) и неявными методами 3-го и 2-го порядка, а также явным методом Эйлера 1-го порядка. Оказалось, что неявные методы в применении к линейным уравнениям (см. рис. 3) существенно менее точно описывают динамику погрешности по сравнению с явными. Поскольку явный метод 1-го порядка дает качественно более близкий (к 4-му порядку) результат, можно сделать вывод, что порядок точности не так сильно влияет на нечеткие вычисления с линеаризацией, как явность метода. В нелинейном же случае (см. рис. 4) не обнаружено различий погрешности решений методами разного порядка, превышающих различие самих решений. Очевидно, более адекватное по сравнению с линейным случаем поведение погрешности связано с применяемым итерационным алгоритмом решения алгебраических уравнений. Также из рис. 4Б видно (особенно при 3-м порядке метода), что разность погрешностей имеет низкочастотные колебания, связанные с набеганием фазы.
3.2. Сопоставление результатов в задаче теории массового обслуживания.
Экономичность метода. Влияние способа описания нечеткости
Ниже метод линеаризации сравнивается, в том числе с точки зрения производительности, с методом многократного решения четких задач на примере нечеткой задачи из области теории массового обслуживания. Эта задача описана в [11,Ошибка! Источник ссылки не найден.], а метод (альтернативный) изложен в [11]. Кроме того, на этом примере результаты расчетов методом линеаризации в гауссовских числах сравниваются с несколькими вариантами интервальных чисел.
Данная задача возникает, например, при земляных работах с одним экскаватором и N грузовиками. После загрузки экскаватором за время 1/ц (ц - скорость обслуживания) грузовики отвозят и разгружают грунт, возвращаются и спустя время 1/у встают снова в очередь к экскаватору. Из-за непредсказуемых изменений параметров ц и V, вероятности возникновения очереди из к грузовиков (обозначаемые ниже через Рк(г), к = 0,1,..^) являются нечеткими. Эти вероятности определяются из следующей системы дифференциальных уравнений:
Р&(г) = - ) + МР1,
Рк (г) = (N - к +1)урк_1(г) - (| + (N - к) V)рк (г) + мрк+1(г), к = 1.^ -1,
& =VРN-1(г) -MРN .
В приведенных ниже расчетах 1/ц принимается равным нечеткому числу с «треугольной» функцией принадлежности с основанием [2;6], а — [9; 11]. В начальный момент времени состояние системы полностью определено: Р^(0) = 1, а для к = 0,1,..N-1 Рк(0) = 0.
На рис. 5 для случая N = 3 сопоставляются результаты решения данной задачи из [11] (на основе сплайн-интерполяции точек, полученных решением четких задач) и результаты метода нечеткой линеаризации, примененного к интервальным числам (которые соответствуют четырем а-уровням разбиения исходного «треугольного» нечеткого числа). Как видно, несмотря на то, что используется стандартная интервальная алгебра (дающая, сама по себе, наиболее широкие интервалы результатов операций), характерная неопределенность решения методом линеаризации меньше неопределенности, вычисленной по минимаксному критерию (но имеет качественно близкую динамику, включая эффект минимума погрешности в окрестности г = 5 с). Меньшая точность расчетов по сравнению с [11] проявляется в симметричности интервалов относительно среднего решения, однако это является свойством использованной интервальной алгебры, а не метода линеаризации.
А. Р2(г) методом линеаризации Б. Р2(г) из [11]
Рис. 5. Сравнение методов при N = 3: вероятность очереди из двух грузовиков Р2(г)
На рис. 6 для случая N = 2 приведена погрешность решения, полученная методом линеаризации с использованием нескольких вариантов интервальной алгебры, гауссовских чисел, а также прямым расчетом 102 четких задач (аналогично [11], только без
интерполяции). На фоне количественного различия между методами (см. положение точки минимума), заметна малая чувствительность результатов к выбору алгебры нечетких чисел.
Рис. 6. Характерная неопределенность вероятности очереди из двух грузовиков p2(t) при N = 2, полученная с помощью разных методов:
• метод сведения к набору четких задач (оранжевая кривая),
• метод нечеткой линеаризации на базе гауссовских чисел (зеленая кривая),
• -//- стандартной интервальной алгебры (синяя кривая),
• -//- интервальных алгебр, уточненных в предположении гауссовской и равномерной функций распределения вероятности (голубая и лиловая кривые)
Для сравнения производительности указанных методов затраты машинного времени измерялись программными средствами, как разность показаний системных часов перед расчетом и по его окончании. Все представленные результаты относятся к персональному компьютеру Pentium III 500 МГц под управлением операционной системы Windows 2000. В связи с оптимизациями операционной системы (в частности, с параллельными операциями), один и тот же расчет как элемент серии занимает меньше времени, чем в отдельности.
При решении задачи обслуживания N = 3 грузовиков методом нечеткой линеаризации с гауссовскими числами затраты машинного времени (за 400 шагов явного метода Рунге-Кутты 4-го порядка) составляют 1.0 с. При использовании различных вариантов интервальной алгебры с представлением нечетких чисел в виде 1 и 4 интервалов эти затраты увеличиваются. В случае стандартной интервальной алгебры затраты составляют 1.5/2.5 с, а в случае ее модификаций (на базе предположений о равномерном и нормальном распределении вероятностей) — 1.8/2.9 с и 1.4/2.1 с, соответственно.
Для оценки аналогичных временных затрат при многократных расчетах четкой задачи предполагается, что каждый нечеткий параметр может быть принят равным одному из M = 5 четких значений (это минимум, необходимый даже для приближенного анализа и даже при использовании сплайн-интерполяции, как это сделано в [11]). Соответственно, два нечетких параметра в силу их независимости дают M2 = 25 комбинаций четких значений, K параметров — MK. Для каждой такой комбинации нужно проводить свой расчет, не говоря уже о дополнительном п-кратном увеличении числа расчетов при аппроксимации функции принадлежности п интервалами. В [11] число п варьировалось от 4 до 8; здесь при оценках затрат оно равнялось 1 (при этом нечеткая задача заменяется на интервальную).
Поскольку в данной задаче число нечетких параметров K = 2, проводился 25-кратный расчет, который в указанном случае занял 3.7 с машинного времени, что в 2.5 раз больше, чем расчет методом линеаризации в стандартных интервальных числах (с п = 1), и в 3.7 раз больше, чем в гауссовских числах.
4. Нечеткое решение алгебраической модели сердца
Ниже коротко рассматривается алгебраическая система уравнений модели сердца (описанной в [14]) и представляются некоторые результаты ее решения в нечетких числах, а также анализ ее чувствительности к исходным данным (который также возможен с помощью метода нечеткой линеаризации). Более подробно эти результаты изложены в [15].
4.1. Параметры модели и ее четкое решение
Рассмотренная в [14] модель сердца имеет 10 параметров (см. таблицу 2), из которых 2 параметра (/ и т^) являются точно измеримыми временными характеристиками работы сердца, 3 безразмерных параметра (К/Е, в, у) описывают индивидуальные и ситуационные особенности эмпирических РК-диаграмм левого желудочка, а 5 параметров (Е, Я1, Я2, АР) характеризуют состояние сосудистых систем и правого желудочка.
Таблица 2. Значения входных параметров модели сердца в норме
Параметр Значение
Частота сердечных сокращений / 65 уд/мин
Длительность систолы т* 0.28 сек
Относительная жесткость желудочка К/Е 0.781 -
Инотропный коэффициент в 1 —
Коэффициент объема сердца У 1 —
Жесткость артерий большого круга Е 0.8 мм.рт.ст/мл
Общее сопротивление большого круга Я1 1.40 мм.рт.ст/(мл/с)
Общее сопротивление малого круга Я2 0.10 мм.рт.ст/(мл/с)
«Сопротивление» правого желудочка Ъ -0.143 мм.рт.ст/(мл/с)
Давление «нежелудочковых насосов» АР 0 мм.рт.ст.
В результате численного решения уравнений модели получались значения 9 переменных, представленных в таблице 3. Формальными переменными решаемой системы были выбраны Р1о„ь Р2ом и Ра, которые отличаются от остальных переменных тем, что при вычислениях правой части уравнения простых итераций используются их значения на предыдущей итерации. Начальное приближение для этих переменных задавалось равным 6, 3 и 82 мм.рт.ст., соответственно; это обеспечивает быструю сходимость (влияние на нечеткое решение средних значений параметров и начальных приближений рассмотрены в
[15]).
Таблица 3. Результирующие переменные модели и их значения в норме
Переменная Значение
Систолический объем Vs 74 мл
Диастолический объем 142 мл
Ударный объем АК 68 мл
Системный кровоток а 73 мл/сек
Систолическое (верхнее) артериальное давление Ps 119 мм.рт.ст
Диастолическое (нижнее) артериальное давление Ра 82 мм.рт.ст
Давление на входе в большой круг Р1т 100 мм.рт.ст
Давление на выходе из большого круга Р1оМ 3 мм.рт.ст
Давление на входе в малый круг Р2т 13 мм.рт.ст
Давление на выходе из малого круга (диастолическое давление в желудочке) Р2оШ 6 мм.рт.ст
4.2. Результаты нечетких расчетов модели и ее чувствительность к параметрам
В рассматриваемых ниже численных экспериментах использовались следующие количественные характеристики нечеткости (неопределенности) параметров модели сердца. При расчетах с полным набором параметров их относительная погрешность бралась равной 15%, с основным набором (параметры К/Е, в, хуже всего известные) — 20%, с одним
нечетким параметром К/Е — 30%, а с одним нечетким сопротивлением большого круга — 40% (физиологический диапазон изменения этого параметра довольно широк). Погрешность вводилась в модель в форме гауссовских чисел, которые обладают простой алгеброй и часто используются при представлении экспериментальных данных (в форме а±с). 25
20 — 15 — 10
0
Нечеткие параметры:
□ Все параметры (15%) ■ К/Е, beta, Rv (20%)
□ К/Е (30%)
□ R1 (40%)
Hi
Т—- --—I-- --—|-- ---1-- --—|-
Рэ Ра Р1 ¡п Р1о1Л Р21п
Рис. 7. Абсолютные погрешности результирующих давлений (мм.рт.ст.)
50%
P2out
40% 30% 20% 10% -Н
0%
Нечеткие параметры:
□ Все параметры (15%) ■ К/Е, beta, Rv (20%)
□ К/Е (30%)
□ R1 (40%)
mllhiH
о.
Q
Ра
P1in
Уб УС! С!У О РБ
Рис. 8. Относительные погрешности переменных модели
На рис. 7 показана гистограмма абсолютных погрешностей (о) результирующих переменных, имеющих одинаковую размерность (давления) при различных наборах нечетких параметров. Более широкое множество переменных представлено на рис. 8 (они имеют разную размерность, поэтому по оси отложены относительные погрешности).
Метод нечеткой линеаризации для решения уравнений с нечеткими параметрами одновременно является инструментом анализа чувствительности этих уравнений к параметрам, хотя и с некоторыми оговорками (чувствительность значения переменной и чувствительность ее погрешности — разные вещи). Ниже этот факт иллюстрируется путем представления коэффициентов су разложения результирующих переменных модели сердца (хг) по исходным данным (¿у): см. формулу (1). Наряду с самими коэффициентами (величинами размерности [переменная]/[параметр]) ниже используются следующие производные понятия, облегчающие сопоставление коэффициентов между собой. Во-первых, это «вклад» параметра в значение переменной sij = са (где ау- - среднее значение ¿у); данная величина имеет размерность переменной. Во-вторых, очень полезным понятием является безразмерный «приведенный вклад» параметра в значение переменной
Ян = su
С
stk\ + \xt о|
V k
\
-1
Z\яь
j
< 1.
Наконец, «нормированный на максимум» вклад Qij больше других подходит для сравнения коэффициентов при одном параметре в разных переменных (для каждой переменной всегда существует один вклад, по модулю равный единице) и не содержит свободного члена
линейной комбинации (из-за больших значений которого приведенные вклады всех параметров в некоторые переменные могут оказаться очень малыми):
Qij = si,
k
-1
: Ы *1-
На рис. 9-11 представлены коэффициенты и вклады параметров в результирующие переменные при расчетах с одним (рис. 9) и со всеми (рис. 10-11) нечеткими параметрами. Жесткость артерий E вычислялась по известной из статистики кривой зависимости жесткости от возраста, поэтому на гистограммах вместо E фигурирует параметр '"age".
На основе представленных результатов можно сделать следующие выводы. Во-первых, величины относительных погрешностей результирующих переменных (см. рис. 8) в большинстве случаев меньше (иногда существенно) относительных погрешностей каждого из нечетких параметров. Особенно важно, что это касается основных результатов модели — артериальных давлений, ударного объема, кровотока: при расчетах с любыми наборами нечетких параметров (характеризующихся погрешностями в 15-40%) эти показатели имеют неопределенность всего 5-10%. И напротив, максимальной неопределенностью обладают наименее значимые в практическом плане результаты: так, относительная погрешность систолического и диастолического объемов (разность которых дает упомянутый выше ударный объем) в некоторых случаях почти достигает 40%, а абсолютная погрешность «низких давлений» (которые также очень сложно измерить на практике) даже превышает погрешность артериальных давлений, имеющих значительно большие четкие значения.
Из рис. 10 видно, что симпатический инотропный эффект (увеличение в) оказывает отрицательное влияние только на низкие давления и диастолический объем (увеличивая остальные, более значимые переменные), а возраст (повышение жесткости артерий и желудочка) существенно увеличивает верхнее артериальное давление и уменьшает значения всех остальных переменных. Рис. 11 иллюстрирует тот факт, что нижнее артериальное давление уменьшается только с возрастом (с ростом остальных параметров увеличивается), а системный кровоток падает также при увеличении системного сопротивления R1 и параметра K/E. Наибольший (положительный) вклад в нижнее артериальное давление вносит R1, в верхнее артериальное давление — K/E, а в системный кровоток — частота сердечных сокращений / 30
20
10
0
-ю -н
-20 -30 + -40 --50 -I--60 -70
Коэффициенты
при параметрах: И R1 ■ beta DRv
Vs Vd dV Q Ps Pa P1in P1 out P2in Рис. 9. Коэффициенты разложения результирующих переменных по параметрам (R1, в, Rv) при расчетах с одним нечетким параметром
0.2 0.15 0.1 0.05 0
-0.05 -0.1 -0.15
Рис. 10. Приведенные вклады некоторых параметров (Р, у, age, f) в результирующие переменные при расчете со всеми нечеткими параметрами
1
0.8
0.6 0.4
0.2
0
-0.2
Рис. 11. Нормированные вклады параметров в основные результирующие переменные Рз, Ра) при расчете со всеми нечеткими параметрами
Все перечисленное полностью согласуется с известными эффектами, в то время как можно привести и неочевидные соотношения между параметрами и переменными. Например, интересным результатом является относительно слабая зависимость основных переменных от неточно известного параметра коэффициенты при нем на порядок меньше коэффициентов при других основных параметрах, а из-за малого четкого значения его вклад меньше других на два порядка (см. рис. 11). Стоит также отметить очень большой вклад частоты пульса в системный кровоток: в расчете со всеми нечеткими параметрами его вклад в 3.5 раза превышает вклад ближайшего «конкурента» по влиянию на кровоток — инотропного коэффициента р.
Вычислительные эксперименты с нечеткой моделью сердца показали экономичность метода нечеткой линеаризации, которая проявляется в алгебраических задачах не в меньшей степени, чем в дифференциальных. Так, 10 (простых) итераций расчета в гауссовских числах с одним нечетким Я1 требует 0.10 с, с основным набором нечетких параметров — 0.15 с, с полным набором — 0.14 с. В то же время, один четкий расчет занимает 0.072 с, а его 5, 125 и 107-кратные повторения, минимально необходимые (см. раздел 3.2) для расчета с указанными наборами параметров, требуют 0.19 с, 1.8 с и 9.5-104 с — соответственно, в 1.9, 12 и 7-105 раз большее время, чем в методе нечеткой линеаризации.
5. Резюме
В работе предложен метод проведения алгебраических операций с зависимыми нечеткими числами на основе хранения линеаризованной истории операций. Метод
Приведенные вклады параметров:
— □ beta ■ aamrna
□ age
— V n_ Df
и i i . LJ . . Щ_ -
и
Vd dV Q Ps Pa PI in P2out
Нормированные вклады в переменные: —
□ G — г>,-
et ш 1 □ Pa
b a gamma £ agí 1 и 5 К/Е f i i Rv L R1
позволяет реализовывать расчеты с нечеткими числами без привязки к конкретным классам задач (алгебраических или дифференциальных) и к вычислительным алгоритмам, предназначенным для решения этих задач в вещественных числах. Кроме того, он позволяет анализировать результаты расчетов на предмет того, какие исходные данные внесли в них наибольший вклад, а также максимально экономичным образом проводить серию однотипных расчетов, в которых нечеткости исходных данных отличаются по величине и форме представления.
В работе представлены результаты расчетов алгебраических и обыкновенных дифференциальных уравнений в интервальных и гауссовских нечетких числах с помощью нескольких стандартных вычислительных алгоритмов в сочетании с методом нечеткой линеаризации. Для дифференциальных уравнений проведено сопоставление полученных погрешностей решения с их оценками через многократное решение соответствующих четких задач, а также с результатами других авторов, пользующихся основанной на этом методикой решения нечетких дифференциальных уравнений. Между ними обнаружено хорошее соответствие, хотя применение неявных численных методов может приводить к качественно неверному поведению нечеткости решения линейных колебательных систем (чего не наблюдается при решении аналогичных нелинейных систем с использованием итерационных методов расчета алгебраических уравнений). Для колебательных уравнений корректно описывается эффект появления низкочастотной составляющей погрешности при нечеткой основной частоте колебаний. Для алгебраической модели сердца показано соответствие нечетких результатов различным физиологическим соображениям, включая результаты анализа чувствительности модели к исходным данным, также проведенного методом линеаризации.
Вычислительные затраты при расчетах методом нечеткой линеаризации (теоретически пропорциональные числу нечетких параметров К) в рассмотренных дифференциальных примерах оказались в несколько раз ниже затрат на многократные расчеты четких задач (теоретически пропорциональное несмотря на малое число К в примерах. С ростом числа параметров относительная экономичность метода линеаризации становится еще более заметной; для модели сердца с 10-ю параметрами она достигает шести порядков.
Литература
1. Алтунин А.Е., Семухин М.В. Модели и алгоритмы принятия решений в нечетких условиях: Монография. - Тюмень: Издательство Тюменского государственного университета, 2000. 352 с. http://www.plink.ru/tnm/cont.htm.
2. Калмыков С.А., Шокин Ю.И., Юлдашев З.Х. Методы интервального анализа. -Новосибирск: Наука, 1986, 222 с.
3. Абрамович Ф.П., Вагенкнехт М.А., Хургин Я.И. Решение нечетких систем линейных алгебраических уравнений LR-типа. // В сб.: Методы и системы принятия решений. -Рига: РПИ, 1987, с. 35-47
4. Захаров А.В., Шокин Ю.И. Алгебраическое интервальное решение систем линейных интервальных уравнений Ах = b и Ах + d = b. // Препринт ВЦ СО АН СССР, N 5, — Красноярск, 1987, 17 с.
5. Семухин М.В. Разрешимость нечетких и интервальных уравнений. Вестник Тюменского государственного университета, вып.2. - Тюмень, ТюмГУ, 1998, с. 23-26.
6. Аленфельд Г., Херцбергер Ю. Введение в интервальные вычисления - М: Мир, 1987, 360 с.
7. Кейн В.М. Оптимизация систем управления по минимаксному критерию. - М: Наука, 1985, 248 с.
8. Базаров М.Б., Шокин Ю.И., Юлдашев З.Х. О построении конечно-разностных интервальных методов для обыкновенных дифференциальных уравнений. // В сб.: Вопросы вычислительной и прикладной математики, ИК АН УзССР, вып.71, 1984, с. 131-144.
9. Маланин В.В., Полосков И.Е. Случайные процессы в нелинейных динамических системах. Аналитические и численные методы исследования. - Ижевск: НИЦ «Регулярная и хаотическая динамика», 2001, 160 с.
10. Schwandt H. A symmetric iterative interval method for system of nonlinear equations. // Computing, 1984, N2, p. 153-164.
11. Oberguggenberger M., Pittschmann S. Differential equations with fuzzy parameters. // Mathematical and Computer Modeling of Dynamical Systems, 5:181-202, 1999.
http://techmath.uibk.ac.at/numbau/publications/98-2.ps
12. Fetz T., Jäger J., Köll D., Krenn G., Lessmann H., Oberguggenberger M., Rieser A., Stark R. F. Fuzzy models in geotechnical engineering and construction management // Computer-Aided Civil and Infrastructure Engineering, 14:93-106, 1999. http://techmath.uibk.ac.at/numbau/publications/97-4.ps
13. Antonelli P., Krivan V. Fuzzy differential inclusions as substitutes for stochastic differential equations in population biology. // Mathematical Biology and Medicine, Vol. 2, World Scientific, Singapore, pp. 165-183., 1996. http://www.entu.cas.cz/krivan/papers/osid92.pdf
14. Евдокимов А.В., Холодов А.С. Квазистационарная пространственно распределенная модель замкнутого кровообращения организма человека // Компьютерные модели и прогресс медицины. М.: Наука, 2001.
15. Евдокимов А. В. Применение метода нечеткой линеаризации к численному моделированию работы сердца // Компьютерные модели и прогресс медицины. М.: Наука, 2003.
