Использование метода недоопределенных вычислений для анализа чувствительности математических моделей Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 4, № 4, 1999

ИСПОЛЬЗОВАНИЕ МЕТОДА НЕДООПРЕДЕЛЕННЫХ ВЫЧИСЛЕНИЙ ДЛЯ АНАЛИЗА ЧУВСТВИТЕЛЬНОСТИ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ

Т. П. КАШЕВАРОВА Российский институт искусственного интеллекта Новосибирский филиал, Россия

In the paper, we consider the use of subdefinite computations technique for the sensitivity analysis of the solutions with respect to variations of parameters in mathematical models. The approach enables one to get guaranteed lower and upper bounds of the values of sensitivity function over finite ranges of the parameters as well as guaranteed bounds of the solution variation. The use of subdefinite computations leads to more accurate estimates of the difference sensitivity function when the solution is not monotonic and unique over the parameter domain.

Введение

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

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

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

Таким образом, методы анализа чувствительности позволяют, с одной стороны, предсказать поведение решения в некоторой окрестности изменения параметров, с другой, —

© Т. П. Кашеварова, 1999.

исходя из требований точности решения модели, оценить допуски в задании эмпирических параметров. Следует отметить, что при использовании функций чувствительности можно эффективно решать задачи управления и оптимизации по параметрам модели. Применение метода недоопределенных вычислений (МНВ) помогает более точно провести анализ чувствительности, так как в этом случае мы получаем гарантированную оценку решения на интервале изменения параметров. Иногда в силу паразитного интервального расширения вычисляемых функций эта оценка получается завышенной, но в данной работе предлагаются способы уточнения интервального решения.

1. Метод недоопределенных вычислений

Метод недоопределенных вычислений разработан в Российском институте искусственного интеллекта в рамках работ по представлению и обработке знаний [1-3]. Рассмотрим применение МНВ для решения системы из т алгебраических уравнений от п переменных:

Р (х) = 0, (1)

где Р = (/*(х), /2(х),... , /т(х)), х = (х\,...,хп) — вектор вещественных переменных, причем некоторые параметры, входящие в функции / (х), могут быть заданы в виде интервалов. Обозначим множество всех решений системы (1) через Р = { х Е К | Р(х) = 0 }, где Р принадлежит параллелепипеду В:

В = { х Е Кп | 1г <хг <гг,1г,гг Е К, г = 1, 2,... ,п }.

МНВ позволяет решить задачу внешнего интервального оценивания, т. е. найти параллелепипед В' такой, что Р С В' С В. Вкратце суть этого метода описана ниже.

Рассмотрим г-е уравнение системы. Так как в общем случае оно содержит не все переменные системы, то его можно записать в виде

/г (хгх , • • • , х1гч ) 0,

где пг < п. Выразив из этого уравнения каждую переменную через другие, получим пг уравнений (называемых функциями интерпретации):

хгх У*! (хг2 , ... , хг„4 ) ,

¿•(г) / \

гз •> 3 У г2 •> ■ ■ ■ 1 "^Ч-п ч ■ ■ ■ ч ) >

хгП1 = (хг1, хг2 , . . . , хггч-1). (2)

Записав в таком виде все уравнения системы (1), получим систему п* +п2+.. .+пт равенств вышеуказанного вида, определяющих в общем случае для г-й переменной системы тг (г < тг < т) функций от разных наборов переменных.

Введем следующие определения и обозначения. Пусть областью значений переменных х1,... , хп является некоторое множество А. Обозначим через *А множество всех подмножеств А и рассмотрим переменные *х, определенные на области *А. Будем называть значения из *А недоопределенными значениями, а переменные *х — недоопределенными переменными. Определим операции над недоопределенными значениями и функции от недоопре-деленных переменных следующим образом. Если 0 — бинарная операция над элементами из А, то ее недоопределенное расширение *0 есть

*а* 0 *Ь := { а 0 Ь | а Е *а, Ь Е *Ь, *а, *Ь Е *А }.

Соответственно недоопределенное расширение */ функции / от п переменных над А есть

7(*Х1,..., *х„) := *{ /(Х1,...,ХП) | х% Е А }.

Будем считать, что все значения переменных системы (1) недоопределены, и недоопреде-ленные значения представим в виде обычных интервалов числовой оси К. Тогда, используя введенные операции и функции для недоопределенных значений, построим следующий итерационный процесс уточнения недоопределенных значений переменных, исходя из представлений (2):

*хк+1 = 7(*хк+1 *хк+1 *хк+1 *хк *хк ) п *к *0 = [1 г ] (3)

<Ау 3 — 33 у •Х'1 ^ •^-'2 ) ' ' ' ) 3_1 ' з+1?***) П'/ 3 ' % — \у%) ' V/

Таким образом, значения переменных в процессе итерационных вычислений представляются последовательностью нерасширяющихся интервалов.

В результате применения итерационного процесса получим либо многомерный параллелепипед *Б, являющийся декартовым произведением *х,к, ] = 1,... ,п, который при = 0 гарантированно содержит все корни системы (1), либо несовместную систему, если = 0. В общем случае не всегда будет самым узким параллелепипедом, содержащим множество решений Р, т. е. можно гарантировать лишь включение Б' С * Д но для некоторых частных случаев [4] доказана сходимость итерационного процесса (3) к единственному решению системы (1).

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

2. Методы анализа чувствительности

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

Г (а, х) = 0, (4)

где а = ( а1, а2, . . . , ак) — вектор параметров, относительно которых изучается чувствительность, х = (х1, х2,..., хп) — вектор решения, Г = ( /1, /2,... , /т), причем в общем случае т = п. Рассмотрим вопросы анализа чувствительности решения х Е Кп в зависимости от параметров а Е Кк. Традиционно этот анализ состоит в вычислении функций чувствительности и матриц чувствительности и в их интерпретации.

Пусть необходимо исследовать чувствительность решения уравнения (4) относительно вариации а,. Функцией чувствительности г-й компоненты решения х» относительно пара-

дх»

метра а,- называют частную производную —— как функцию параметра а. Соответственно

да

уравнения для определения функций чувствительности получают дифференцированием (4) по а, :

дЛ + М дЖ = 0 I = 12 П (5)

да, дх» да, ' >>•••>•

джг

Решая систему алгебраических уравнений (5) относительно —— совместно с системой

да,

уравнений (4) в точке а = а0, получим значение функции чувствительности в этой точке.

Если рассматривать только малые отклонения а, от а0, то для оценки хДа° + 8) можно использовать формулу Тейлора, ограничиваясь первыми членами разложения:

к г\

хДа° + 8) и хг(а°) +^ (а0) 8, '

3 = 1 3

где 8з > 0, 8 = (, 82 , ■ ■ ■ , 8к). В этом случае элементы матрицы чувствительности первого порядка определяются следующим образом:

Ь = дж

3 = да, ■

Когда рассматривается изменение х» в зависимости от а,, все остальные параметры а^, I = , фиксируются. Для каждого параметра а° требуется решать систему 2т уравнений. Матрица чувствительности может быть использована для получения решения линеаризованного уравнения, определяемого как

х(а) и х° + В( а - а°),

хотя проводимый таким образом анализ модели будет наиболее точным только при малых значениях 8 и приведет к значительным ошибкам, когда имеет место смена монотонности решения х» на интервале [а° , а° + 8,].

дхг 1 дхг 2

К'-аЧ |< дО" = ^ и ШаХ1а^'даТ = ^'

Из (4)-(5) можно получить оценки для min^.-„о|<¿ = ßj и max|„.-o0|<¿ —- = ß

3

так как по ним выписывается гарантированная оценка для решения на всем интервале

[а° , а° + 8, ] :

х° + (а, — а°) < ж» < х° + (а, — а°) , Уа, € [а° , а° + 8,]■

В некоторых случаях поиск определенной выше функции чувствительности становится затруднительным. Когда не требуется высокой точности локального исследования решения в окрестности точки а°, можно ограничиться разностной оценкой функции чувствительности, т.е. воспользоваться простым и весьма популярным в вычислительной практике подходом [5]. При этом точная производная заменяется разностным отношением. Именно, вычисляется решение х°, г = 1 , ■ ■ ■ , п, такое, что

^( а° , а° , ■ ■ ■, а° , ■■■, ак , х° , ■ ■ ■, хЩ) = 0 , (6)

затем а3- изменяется на 83 и находится новое решение х3 из системы уравнений

^(а0,а2,... ,аз-1 ,а0 + 83,а°+1,..., а^, х 1,..., х^) = 0. (7)

Разностным аналогом функции чувствительности является следующая величина:

х'3 _ <х0 •

п •X' п {-г •Л-' 1

83 да3

(которая зависит также от 8). Разностная матрица чувствительности В = (63) при таком подходе строится следующим образом:

3 о

63 = , (8)

83

где х0 — компонента решения векторного уравнения (6), х3 — решение уравнения (7),

дхг

% = 1,... , п, ^ = 1,... , к. Если 8 ^ 0, то получим частные производные 63 = —— в точке

аа3

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

Можно ввести обобщенную характеристику чувствительности модели относительно параметра а3, используя некоторую норму 11 • 11 в пространстве Кга и полагая О = || (х_х0)/8311. Например, если взять в качестве || • || евклидову или чебышевскую норму, то получим

О? = шах |. (9)

Евклидова норма дает некоторое усредненное значение чувствительности всего вектора решения х к параметру а3, норма (9) позволяет выделить ту компоненту решения, на которой достигается максимальное изменение решения. В дальнейшем ограничимся рассмотрением нормы (9). При этом имеет смысл запоминать и номер %, на котором достигнут максимум. В этом случае можно ответить на вопрос, какая компонента решения наиболее чувствительна и к какому параметру, если исследуемых параметров несколько.

Преимущество использования МНВ заключается в том, что он позволяет получить гарантированные нижнюю (хт;п) и верхнюю (хтах) оценки для решения вектора х на интервале [а0, а0+83]. Для этого необходимо вместо (6), (7) рассмотреть одно интервальное уравнение

^(а?, а2,..., [а0, а0 + 83],..., а£,х?,..., хП) = 0. (10)

Обозначим [х3 т;п,х3 тах] — компоненты вектора решения, полученные методом недо-определенных вычислений из (10); [х3 т;п,х3 тах] — точные нижняя и верхняя границы компонент решения (10) на интервале [а0, а0 + 83]. Так как МНВ находит интервал, гарантированно содержащий все решения системы, то для этих двух решений будем иметь

[х>3 х'3 ] ( [х^3 <х3 ]

г тт> г тах] — и^г тт> г тах]"

Если [xj min, xj max] намного уже интервала [Xj min, Xj max], то последний можно также сузить, исключая интервалы, не содержащие решения. Для этого МНВ применяется к (10) с дополнительным ограничением

<Xj - [xX^ I ]

•^i b^i min> ^i max TU1)

где dj > 0, и если получим, что система несовместна, то нижняя граница сдвигается на величину dj. Для совместной системы интервал [a0,a0 + ój] сузится, и мы получим окрестность ají, в которой находится уточненная нижняя граница xj. Уменьшая dj, можно добиться желаемой точности. Аналогичная процедура совершается для правой границы с той лишь разницей, что в качестве дополнительного ограничения выбирается

<Xj - ['Xj _ Гj XXj ] Г j >> 0

i i min i max

Таким образом, при нахождении xj min, xj max определяются и соответствующие значения

aji,aj-2, для которых достигаются эти экстремальные значения. При этом величина

jj

Ü2 — max ; i max i min

i

ój

будет оценкой сверху для величины Dí, а

D3 — max

jj x_x

i max i min

i i aj2 _ ají

будет оценкой сверху для Д2. Следовательно, получаем цепочку неравенств:

Д < В2 < Дз.

Если рассмотреть величину

jj i max i min

D4 — max

i

aj2 _ ají

то она будет оценивать Д сверху, т.е. Д < Д4, однако ничего нельзя сказать о том, как соотносятся Д и Д4.

Матрица Ф = (ф^) такая, что

jj x_x

фг] = -^^, (11)

aj2 — ajl

будет более точно по сравнению с (8) отражать максимальную чувствительность решения к изменению параметров на интервале [а0, а0 + 8]. Таким образом, использование МНВ для анализа чувствительности позволяет более точно оценивать значение функции чувствительности для конечных значений 8.

3. Пример анализа чувствительности модели народонаселения

В качестве примера проведем анализ чувствительности модели народонаселения, разработанной в институте ЕАШ (Германия) [6]. Эта модель позволяет прогнозировать численность населения и оценивать стоимость реализации схемы его пенсионного обеспечения. Предполагается, что все население страны разбито на две социальные группы: в первой (группе ) введено гарантированное пенсионное обеспечение в старости, вторая (группа ) такой поддержки не имеет. Переход населения из группы А в группу В происходит в соответствии со значением некоторого параметра — уровня охвата пенсионной схемой. Для простоты рассмотрение ограничивается только женской частью населения.

Математическая модель, описывающая данный процесс, выглядит следующим обра-

зом:

16;

AU,-k = Aj-itk(1 - AMj-i)5, j = 2, AUi>fc = 0;

AUi7,fc = Ai6,fc(1 - AMi6)5 + Ai7>fc(1 - AM17)5;

10

A

i,fc+i

AFj

- AM,

1n(1 - AM,)

) (5 - 10AM, + 10AM,2 - 5AM,3 + AM,4);

,=i ----v- 3J

Aj,k+i = AUjk(1 - acceptk * distrib,), j = 2,..., 17; i0

Bi,k+i = BGk = 2 ^^ Bj,kBFj

- BMj

1n(1 - BM,),

-) (5 - 10BM, + 10BM,2 - 5BM,3 + BMj1), )

j=i ----v j

Bj,k+i = Bj-ik(1 - BMj-i)5 + AUjk ■ acceptk ■ distrib,, j = 2,... , 16; Bi7,fc+i = Bi6>fc(1 - BMi6)5 + Bi7>fc(1 - BMi7)5 + AUi7,fc ■ acceptfc ■ distribi7;

i7

pAfc+i=^ a.

j,k+i;

j=i i7

pBk+i = Y, B.

j,k+i;

j=i

i7

Pfc+i = Pk + 2 ■ p ■ J] Bj;fc+i(

j=i4

BM,

log(1 - BM,)

)(5 - 10BM, + 10BM,2 - 5BM,3 + BM,4);

acceptk = 1 - (1 - a)

pi = 0; 1

.k - b.

1 + exp(-)

c

k = 1, 2, 3,... .

Исследуем чувствительность переменных accept (характеристики уровня рождаемости) , pA (численности населения в группе A), pB (численности населения в группе B), P (суммарной стоимости пенсионного обеспечения) к изменению параметров a, b, c,p в точке ( 0.1, 10, 1.5, 100). Для вектора 6 выберем значение 6= (0.05, 0.5, 0.1, 5). Построим матрицы чувствительности по формуле (11) в выбранной точке для значений времени t = 50 и t = 200 лет (вектор решения (accept,pA,pB, P), вектор параметров (a,b, c,p)):

1

2

t = 50 лет

( 0.33924 -0.141362 -0.08491

-230.5 106.54 V 1537.56

t = 200 лет

(

0 0

11.462 -7.869 -1.94

0 0

-16.64 9.163 15.33

0 0

0 0 0

1.5758

0 0

\

-216.54 20.4 -1.38 0 V -1156600 80040 -34300 11030 )

Анализируя матрицу чувствительности для периода t = 50 лет, можно сделать следующие выводы.

1. Самое существенное влияние на абсолютно все исследуемые компоненты решения оказывает параметр а. Именно его малые изменения приводят к сильному изменению всех выбранных нами для анализа компонент решения. Поэтому к выбору значений параметра а необходимо относиться с особой тщательностью.

2. Параметр p не оказывает абсолютно никакого влияния на переменные accept, pA,pB. Матрица чувствительности, построенная для периода t = 200, существенно корректирует выводы, сделанные при t = 50. Так, параметр а остается ведущим для pB,P, но компоненты решения accept, pA выходят на стационарный режим начиная с некоторого момента времени. Чтобы точнее определить это значение времени, нужно провести серию расчетов, дробя интервал t = [50, 200]. Не приводя полностью всех коэффициентов матрицы чувствительности, укажем, что через 100 лет значение accept выходит на стационарное значение и не зависит от параметров а и b.

Используя метод недоопределенных вычислений, можно получать гарантированные оценки решения для изменения параметров модели в некотором интервале. Приведем интервальное решение для значений параметров а = [0.1, 0.15], b = 10, c = 1.5, p = 100 для времени t = 50 лет:

accept = [0.694681, 0.711643]; pA = [109.374,120.899]; pB = [36.31, 46.35]; P = [157.567, 234.445].

4. Заключение

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

Список литературы

[1] НАРИНЪЯНИ А. С. Недоопределенность в системах представления и обработки знаний. Изв. АН СССР. Техн. кибернетика, №5, 1986, 3-28.

[2] BABICHEV A. B., KADYROVA O.B., KASHEVAROVA T.P., LESHCHENKO A. S., SEMENOV A. L. UniCalc, a novel approach to solving systems of algebraic equations. Int. Comput., No. 2, 1993, 29-47.

[3] SEMENOV A., KASHEVAROVA T., LESHCHENKO A., PETUNIN D. Combining various techniques with the algorithm of subdefinite calculations. In "Proc. of PACT'97". London, 1997, 287-305.

[4] КАШЕВАРОВА Т.П., СЕМЕНОВ А. Л. Некоторые вопросы сходимости метода недо-определенных вычислений. В "Проблемы представления и обработки не полностью определенных знаний". Москва — Новосибирск, 1996, 31-37.

[5] ОРАН Э., Борис Дж. Численное моделирования реагирующих потоков. Мир, М., 1990.

[6] STROHM J. Population Model, Draft. FAW, May 1997.

Поступила в редакцию 2 ноября 1998 г.