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

Данный метод применяется для статистического анализа нелинейных динамических систем, математическое описание которых может быть представлено в виде [6] :

(6.1)

,

где- - n - мерный вектор неслучайных начальных условий; - к - мерный вектор центрированных некоррелированных между собой случайных величин, для которых известны моментные характеристики до q - порядка включительно,

т.е. .

Требуется найти моментные характеристики выходных величин вида:

В основе метода лежит гипотеза об известном характере зависимости между входными случайными параметрами и выходными величинами. Для рассматриваемого случая - это возможность ее представления в виде ряда Тейлора по случайным параметрам V₁,: V_m до q - порядка включительно.

Пусть решения системы уравнений (6.1) являются функциями времени t и случайных величин V_1,:,V_m :

(6.2)

Тогда р -я степень решения (6.2) может быть представлена в виде:

, (6.3)

где

Разложим функцию в ряд Тейлора по случайным параметрам V_1,:,V_m в окрестности их нулевых математических ожиданий, ограничиваясь членами q -й степени :

(6.4)

где - значения функции при нулевых значениях случайных величин, - частные производные от функции по случайным величинам, вычисленные при их нулевых значениях. Применив к (6.4) операцию отыскания математического ожидания, получим искомую вероятностную характеристику

. (6.5)

Однако решение (6.2) обычно неизвестно, поэтому найти частные производные в (6.5) затруднительно. Для устранения этой проблемы выбирается N различных комбинаций значений для входных случайных величин, называемых в дальнейшем эквивалентными возмущениями :

Решив N раз систему (6.1), получим N решений . После подставления эквивалентных возмущений и соответствующих им решений в (6.4), умножения обеих частей полученных N равенств на пока неизвестные коэффициенты и суммирования полученных результатов записывается следующее соотношение :

.(6.6)

Сравнивая выражения (6.5) и (6.6), можно сделать вывод, что сумма Q будет равна , если выполнены условия

(6.7)

где .

В результате решения системы уравнений (6.7) находятся неизвестные .

Исходная система уравнений (6.1) интегрируется N раз для найденных комбинаций эквивалентных возмущений, а затем рассчитываются и находятся искомые вероятностные характеристики по соотношениям вида

(6.8)

Найдем решение системы (6.7) для случая q=2 (соотношения метода эквивалентных возмущений для q = 3 и q = 5 приведены в [6]).

Пусть для системы случайных величин V₁,:,V_m известны вероятностные характеристики

,

(6.9)

Тогда система (6.7) примет вид

(6.10)

Положим N=m+2 и выберем значения эквивалентных возмущений x _rS согласно нижеприведенной табл. 1, добиваясь обращения в нуль возможно большего числа коэффициентов a _S.

Для выбранного числа комбинаций и принятой матрицы эквивалентных возмущений из решения системы (6.10) следует

(6.11)

Таким образом, главные вероятностные характеристики выхода динамической системы M[x_i] и D[x_i] могут быть найдены по соотношениям :

(6.12)

(6.13)

Аналогичные формулы можно получить и при числе членов разложения в ряд Тейлора q > 2.

В заключение приведем этапы статистического анализа нелинейной динамической системы с помощью метода эквивалентных возмущений :

1).Выбор значений эквивалентных возмущений согласно табл. 1.

2).Численное решение системы дифференциальных уравнений (6.1) (m+2) раза для комбинаций значений случайных параметров, равных соответствующим эквивалентным возмущениям.

3).Вычисление искомых вероятностных характеристик по формулам типа (6.12), (6.13).

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

Пусть решение системы (6.1) имеет вид

(6.14)

Требуется найти оценку математического ожидания некоторой функции от фазовых координат системы в момент времени t (под эту формулировку подходят любые вероятностные характеристики) :

(6.15)

Так как явный вид выражения (6.14) неизвестен, то можно найти значения решений (6.14) путем численного интегрирования (6.1) на ЭВМ при определенных значениях случайных параметров V₁:V_m, которые называются узлами интерполирования.

Пусть для первой случайной величины V₁ выбраны определенным образом q₁ вариантов узлов интерполирования ; для произвольной j -й случайной величины V_j выбраны q_j вариантов узлов интерполирования ; для последней m -й случайной величины V_m выбраны q_m вариантов узлов интерполирования .

Тогда общее число расчетных вариантов .

Представим функцию Ф(...) приближенно с помощью интерполяционного полинома Лагранжа :

(6.16)

где - многочлен степени q_j относительно переменной V_j; - производная от по переменной V_j, вычисленная в точке .

Тогда искомую характеристику (6.15) М[Ф(...)] можно представить с учетом (6.16), поменяв местами порядок интегрирования и суммирования, в виде

. (6.17)

Здесь

(6.18)

называются числами Кристоффеля.

В [7] доказано, что если в качестве узлов интерполирования выбрать корни ортогональных полиномов по весу, равному плотности распределения случайной величины p(V), то при использовании n узлов интерполирования данный метод дает точные значения в классе многочленов всех степеней до степени q=2n-1 включительно. При этом для каждой плотности распределения вероятностей, являющейся весовой функцией, существует единственная система ортогональных многочленов [7].

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

Пусть случайная величина V распределена равномерно на интервале [a,b]. Тогда с помощью линейного преобразования вида

(6.19)

случайную величину V можно получить из случайной величины l , распределенной равномерно на интервале [-1, 1] с плотностью .

. (6.20)

Числа Кристоффеля для узлов, выбранных в точках, соответствующих корням многочленов Лежандра степени (n+1), определяются по формуле

. (6.21)

Значения корней l _k ортогональных многочленов Лежандра до 8 порядка включительно (q<=8) для случая закона равномерного распределения плотности вероятности на интервале [-1, 1] приведены в табл. 2. Их следует пересчитывать для конкретных параметров равномерного закона распределения случайной величины по (6.19).

. (6.22)

Числа Кристоффеля r _к рассчитываются заранее и для q<=8 также приведены в табл. 2, при этом они выбираются из таблицы без пересчета.

Пусть случайная величина V распределена по нормальному закону с параметрами M_V, D_V. Тогда с помощью линейного преобразования вида

(6.23)

случайную величину V можно получить из нормально распределенной случайной величины l с Мl = 0, Dl = 1.

Весовой функции вида соответствует система ортогональных многочленов Эрмита на интервале (-¥ , +¥ ), которые можно вычислить по рекуррентным формулам [7] :

(6.24)

Значения корней ортогональых многочленов Эрмита l _к для веса р_Г(l ) приведены в табл. 3 и пересчитываются для случайной величины V по (6.23):

. (6.25)

Числа Кристоффеля r _к не зависят от параметров нормального закона распределения и выбираются также из табл. 3 без пересчета. Таблицы, подобные табл. 2 и 3, могут быть рассчитаны и для других законов распределения случайных величин, например, для экспоненциального закона приведены в [7].

для равномерного распределения вероятностей

для нормального распределения вероятностей.

Опишем этапы статистического анализа динамической системы при использовании интерполяционного метода:

1.Для заданных порядков интерполяционных многочленов q₁,:,q_m по случайным величинам V₁,:,V_m соответственно находятся по таблицам (типа таблиц 2, 3) узлы интерполирования и числа Кристоффеля .

2.По формулам пересчета типа (6.19), (6.25) рассчитываются значения узлов интерполирования для каждой случайной величины - , а также числа Кристоффеля вида (6.18).

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

4.Вычисляется искомая вероятностная характеристика по соотношению (6.17).

Трудоемкость метода можно существенно сократить, если при формировании расчетных вариантов контролировать значение выражения .

Если начинает выполняться неравенство

< d , (6.26)

где d - заданная допустимая погрешность вычислений, то выборка узлов интерполирования прекращается и в дальнейших вычислениях участвуют лишь выбранные узлы [7].

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

, (6.27)

, (6.28)

(6.29) .(6.30)

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

Метод статистических испытаний (или метод статистического моделирования, или метод Монте-Карло) является наиболее распространенным на практике методом статистического анализа нелинейных динамических систем. Его достоинствами являются наглядная вероятностная трактовка, способствующая быстрому практическому усвоению; универсальность применения к исследованию точности любых динамических систем; простая вычислительная схема; существенно упрощающая программирование; устойчивость результата по отношению к возможным ошибкам при проведении отдельных моделирований; простая оценка точности получаемых результатов.

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

, (6.31)

. (6.32)

Далее многократно повторяются реализации алгоритма (выполняется моделирование динамической системы, при которых случайные факторы V, x , x₀ заменяются их реализациями, полученными случайным образом) и формируется статистическая выборка (второй этап) :

, (6.33)

, (6.34)

где , , - реализации случайных факторов V, x , x₀ соответственно, N - объем выборки реализаций.

Затем обрабатывается полученная статистическая выборка реализаций методами математической статистики и находятся искомые вероятностные характеристики выхода динамической системы (третий этап) :

. (6.35)

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

Так как в (6.35) оценка является суммой случайных величин, равнозначных в вероятностном отношении, то согласно центральной предельной теореме Ляпунова закон распределения случайной величины будет при N --> ¥ приближаться сколь угодно близко к нормальному закону распределения. Поэтому в качестве меры относительной погрешности оценки можно принять величину [7]

(6.36)

Тогда можно записать следующую вероятностную оценку погрешности [7] :

(6.37)

где e - заданная верхняя граница погрешности d _N ; - вероятность не превышения погрешностью d _N фиксированной границы e ; Ф(...)- функция Лапласа.

Из (6.37) следует, что если гарантийную вероятность принять равной Р_г = 0,99, то формула для определения необходимого числа опытов N для обеспечения требуемой точности e с вероятностью Р_г запишется в виде

,

откуда

. (6.38)

Таким образом, определение оценки вероятностных характеристик по методу статистических испытаний требует проведения большого объема моделирований (например, достижение 5% -й точности требует проведения не менее 400 испытаний), что является главным недостатком упомянутого метода.

Для снижения трудоемкости метода статистических испытаний разработаны различные его модификации. Одна из наиболее простых и в то же время эффективных модификаций основана на дополнительной обработке выборки входных случайных факторов и случайных выходных контролируемых характеристик [11].

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

Если математическая модель динамической системы преобразована к виду (6.1), то уточненная оценка вероятностной характеристики M*[y_i] определяется по соотношению

. (6.39)

Здесь M[V] - точное значение вектора математических ожиданий входных случайных величин V (известно по условию задачи);

- вектор оценок математических ожиданий входных случайных величин V, найденный в результате статистической обработки выборки реализаций ;

- вектор оценок выходных контролируемых характеристик, найденный в рамках метода статистических испытаний;

- погрешность в определении оценки вероятностных характеристик для входных случайных факторов по методу статистических испытаний;

K_VV - автоковариационная матрица для векторной случайной величины V, которая в (6.39) заменяется своей оценкой K_VV(N), найденной в результате обработки выборки входных случайных факторов

;

- матрица смешанных моментов между y_i и V, отражающая их статистическую взаимосвязь, которая в (6.39) заменяется своей оценкой , найденной в результате обработки выборки входных случайных факторов и выходных контролируемых характеристик :

.

Использование оценки (6.39) позволяет получить более высокую точность при том же числе испытаний по сравнению с методом статистических испытаний или сократить трудоемкость при сохранении погрешности результата.

Оценка снижения трудоемкости (или повышения точности результата) по i -й контролируемой характеристике может быть найдена по формуле

, (6.40)

где - квадрат множественного коэффициента корреляции между y_i и V; - оценка дисперсии контролируемой характеристики y_i.

Из (6.40) следует, что если =0 (что соответствует отсутствию статистической взаимосвязи между входными случайными параметрами и выходными контролируемыми характеристиками), то выигрыш от использования уточненной оценки не наблюдается.

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