ОБРАБОТКА ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ С ИСПОЛЬЗОВАНИЕМ СРЕДЫ SIMINTECH

Необходимость обработки и анализа большого объёма экспериментальных данных привело к появлению разнообразных по функциональному наполнению, удобству в использовании и уровню визуализации статистических пакетов, позволяющих выявить явные и скрытые закономерности. Наиболее известными среди них являются Statistica, EViews и SPSS. С одной стороны, применение статистических пакетов значительно упрощает и ускоряет процесс обработки, с другой – не даёт пользователю в полной мере самостоятельно реализовывать и корректировать алгоритмы для полного понимания вычислительного процесса.

Данная работа демонстрирует подход к решению задач обработки экспериментальных данных с использованием среды динамического моделирования технических систем SimInTech. Цель использования данной среды заключена в переходе от использования статистических пакетов для решения подобного типа задач. Далее рассмотрим последовательность действий при реализации парного сравнения по степени автономности полёта 10 различных образцов управляемых авиационных бомб (УАБ), которую иначе можно сформулировать как задачу анализа различий между двумя независимыми выборками.

Основная проблема заключается в том, что степень автономности УАБ является случайной величиной, поскольку зависит от таких факторов, как: условия сброса, внешние ветровые возмущения, аэродинамические характеристики УАБ. Поэтому необходимо проводить сравнение выборок на основе некоторой числовой характеристики.

Подобного рода задачи решаются одним из методов проверки статистических гипотез: параметрическим методом, то есть необходимо доказать, что распределение набора значений степени автономности сравниваемых УАБ подчиняется нормальному закону с одинаковыми значениями генеральных дисперсий.

Степень автономности полёта является одним из тактических показателей, используемых для оценки технического уровня образцов УАБ, и определяется после сброса бомбы по формуле:

                                                                                                                (1)

где: T_{УПРиНАВ} – суммарное время управления и наведения УАБ на цель;

T_Б – общее время полета УАБ после сброса с авиационного носителя.

В ходе испытаний для каждого из 10 образцов УАБ было произведено 50 бомбометаний по точечной наземной цели и получены значения степени автономности, приведенные в таблице 1:

Таблица 1.

Значения степени автономности полёта

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

Для начала необходимо определить экспериментальную вероятность попадания в i-й интервал, где i – порядковый номер интервала. Количество интервалов варьируется и выбирается самостоятельно. В рамках рассматриваемого примера был выбрано 10 равных по величине непересекающихся интервалов.

Для вычислений указанной выше вероятности в SimInTech реализована схема автоматики «Определение экспериментальной вероятности попадания в i-й интервал» (рисунок 1).

Рисунок 1. Схема «Определение экспериментальной вероятности попадания в i-интервал»

Блок «График Y от X» [1, с. 405] в данной схеме представляет собой гистограмму распределения, отражающую графическую зависимость вероятности попадания элементов выборки в соответствующий интервал группировки. Входным сигналом для первого блока «Язык программирования» [1, c. 381] является блок «Линейный источник» [1, c. 376], представляющий собой функцию y(t) = a+b·t, где a = 0, b = 1/h, h – шаг интегрирования в секундах, t – текущее модельное время в секундах.

Выполнение расчётов осуществляется при помощи скрипта проекта, который содержит следующий код:

h = 1e-3; // шаг интегрирования

X = [0.63, 0.47, 0.39, 0.44, 0.74, 0.36, 0.40, 0.49, 0.41, 0.41, 0.52, 0.67, 0.45, 0.34, 0.34, 0.50, 0.49, 0.43, 0.44, 0.55, 0.27, 0.59, 0.43, 0.54, 0.42, 0.36, 0.52, 0.60, 0.51, 0.52, 0.54, 0.35, 0.71, 0.77, 0.69, 0.52, 0.59, 0.42, 0.61, 0.51, 0.31, 0.46, 0.51, 0.49, 0.52, 0.58, 0.45, 0.57, 0.51, 0.65];

n = cols(X); // количество элементов массива

n_int = 10; // количество интервалов

int = (max(X)-min(X))/n_int; // величина интервала

t_end = h*(n-1); // время окончания расчёта

Расчёт этой и всех последующих схем выполняется с постоянным шагом h при помощи метода Рунге-Кутта 4-5 порядка.

Значения массива X задаются в зависимости от исследуемого образца УАБ из таблицы 1.

Первый блок «Язык программирования» отвечает за извлечение значения i-го элемента из исходного массива X с помощью встроенной функции extract [2]:

input i;

xi = extract(X, [1], [i]);

output xi;

Второй подобный блок содержит код и предназначен для подсчёта количества попаданий элементов выборки в каждый из i-х интервалов:

initialization

var ni = (n_int)#0;

  var i: integer;

end;

input xi;

if goodstep and (xi > 0) and (xi <= min(X)+int) then ni[1] = extract(ni, [1])+1;

for (i=2,n_int)

if goodstep and (xi > max(X)-(n_int-i+1)*int) and (xi <= max(X)-(n_int-i)*int) then ni[i] = extract(ni, [i])+1;

output ni;

Запись var ni = (n_int)#0 соответствует созданию «нулевого» массива с количеством элементов равным количеству интервалов, goodstep – флаг «хорошего» шага (тип данных – boolean).

Выход данного блока делится на количество элементов массива X, что соответствует формуле:

                                                                                                                          (2)

где: n_i – число попаданий значений элементов выборки в i-й интервал;

N – количество элементов в выборке.

Значения математического ожидания m_x и среднеквадратического отклонения σ_х получены в результате моделирования данной схемы посредством встроенных функций mean(X) и rms(X). Далее они будут использованы при определении значений функции Лапласа.

В результате моделирования представленной схемы получены значения экспериментальные вероятности, которые для каждого образца УАБ приведены в таблице 2.

Таблица 2.

Экспериментальные вероятности образцов УАБ

В качестве результата моделирования схемы (рисунок 1) представлена гистограмма распределения для первого образца УАБ (рисунок 2).

Рисунок 2. Гистограмма распределения

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

                                                                                                          (3)

где:  для i = 2..(n_int–1) (i – индекс границы интервала, n_int – количество интервалов),  для i = 1,  для i = n_int.

Схема «Определение значения функции Лапласа» представлена на рисунке 3. Конечное время расчёта равно 100 секундам.

Рисунок 3. Схема «Определение значения функции Лапласа»

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

h = 1e-5; // шаг интегрирования

X = [0.63, 0.47, 0.39, 0.44, 0.74, 0.36, 0.40, 0.49, 0.41, 0.41, 0.52, 0.67, 0.45, 0.34, 0.34, 0.50, 0.49, 0.43, 0.44, 0.55, 0.27, 0.59, 0.43, 0.54, 0.42, 0.36, 0.52, 0.60, 0.51, 0.52, 0.54, 0.35, 0.71, 0.77, 0.69, 0.52, 0.59, 0.42, 0.61, 0.51, 0.31, 0.46, 0.51, 0.49, 0.52, 0.58, 0.45, 0.57, 0.51, 0.65];

mx = 0.4998; // математическое ожидание элементов массива

sigmax = 0.112; // среднеквадратичное отклонение элементов массива

n_int = 10; // количество интервалов

int = (max(X)-min(X))/n_int; // величина интервала

i = 1; // индекс границы интервала (от 1..n_int+1)

if (i >= 2) and (i <= n_int) then ti = (min(X)+(i-1)*int-mx)/sigmax else

if i < 2 then ti = -mx/sigmax else ti = 1/sigmax;

if ti > 0 then t_end = ti else t_end = abs(ti);

Циклы c логическим if соответствуют определению верхнего предела интеграла из формулы (3).

Первый блок «Язык программирования» останавливает расчёт схемы по достижению верхнего предела интеграла, во втором подобном блоке определяется знак функции Лапласа.

Таким образом, изменяя индекс границы интервала i, были получены значения функции Лапласа для каждого образца УАБ.

Таблица 3.

Значения функции Лапласа для образцов УАБ

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

Как упоминалось ранее, для того чтобы рассматривать сравниваемые образцы УАБ как независимые выборки, необходимо доказать, что распределение значений в каждом из сравниваемых образцов является нормальным.

Для данной проверки была построена схема «Проверка выборки на нормальность распределения (критерий χ-квадрат)», приведенная на рисунке 4:

Рисунок 4. Схема «Проверка выборки на нормальность распределения (критерий χ-квадрат)»

Скрипт проекта содержит следующий код:

h = 1e-3; // шаг интегрирования

X = [0.63, 0.47, 0.39, 0.44, 0.74, 0.36, 0.40, 0.49, 0.41, 0.41, 0.52, 0.67, 0.45, 0.34, 0.34, 0.50, 0.49, 0.43, 0.44, 0.55, 0.27, 0.59, 0.43, 0.54, 0.42, 0.36, 0.52, 0.60, 0.51, 0.52, 0.54, 0.35, 0.71, 0.77, 0.69, 0.52, 0.59, 0.42, 0.61, 0.51, 0.31, 0.46, 0.51, 0.49, 0.52, 0.58, 0.45, 0.57, 0.51, 0.65];

n_int = 10; // количество интервалов

n = cols(X); // количество элементов массива

pi_exp = [0.04, 0.1, 0.12, 0.16, 0.26, 0.08, 0.1, 0.06, 0.04, 0.04]; // экспериментальное значение вероятностей

F = [-0.5, -0.4458, -0.3768, -0.2619, -0.1049, 0.0716, 0.2346, 0.3584, 0.4357, 0.4754, 0.5]; // значение функции Лапласа

chi2_theor = 14.0671; // квантиль распределения хи-квадрат (Pd = 0.95, k = n_int-3 = 7 - число степеней свободы)

t_end = h*(n_int-1); // время окончания расчёта

здесь вектор F – вектор значений функции Лапласа для соответствующего образца УАБ из таблицы 3.

Проверка проводится согласно критерию согласия χ-квадрат, расчёт которого производится в блоке «Язык программирования» согласно формуле [3, с. 67]:

                                                                                                       (4)

где:  – найденное экспериментально значение критерия χ-квадрат;

N – количество элементов в выборке;

 – теоретические значение i-й вероятности;

 – экспериментальные значение i-й вероятности, приведённые в таблице 2.

Данная формула реализована в блоке «Язык программирования» при помощи нижеследующего кода:

initialization

  var pi_theor = (n_int)#0; // теоретическое значение вероятностей

  var s = (n_int)#0; // компоненты суммы для определения экспериментального значения критерия хи-квадрат

end;

input i;

pi_theor[i] = extract(F, [i+1])-extract(F, [i]);

s[i] = (((extract(pi_exp, [i])-extract(pi_theor, [i]))^2)/(extract(pi_theor, [i])));

chi2 = n*sum(s); // экспериментальное значение критерия хи-квадрат

output chi2;

Выборка значений исходного массива Х будет принадлежать нормальному распределению лишь в том случае, если экспериментальное значение критерия χ-квадрат не превышает теоретического значения критерия взятого при значении доверительной вероятности P_д = 0.95 со степенями свободы n_int–3, то есть .

В таблице 4 приведены экспериментальные значения критерия χ-квадрат для каждого из 10 образцов УАБ.

Таблица 4.

Экспериментальные значения критерия χ-квадрат

Проанализировав данную таблицу, можно увидеть, что не все элементы выборки каждого из образцов УАБ подчиняются нормальному распределению. В рамках дальнейшего анализа будут рассмотрены только образцы с номерами 1, 3, 7, 8, 9, 10, которые удовлетворяют условию нормальности распределения.

Далее для отобранных образцов необходимо проверить гипотезу о равенстве дисперсий. Для этого реализована схема «Проверка гипотезы о равенстве дисперсий с использованием критерия Фишера», приведенная на рисунке 5.

Рисунок 5. Схема «Проверка гипотезы о равенстве дисперсий с использованием критерия Фишера»

Для проверки гипотезы о равенстве дисперсий используется критерий Фишера, который рассчитывается в блоке «Язык программирования» в соответствии со следующей формулой [3, с. 74]:

                                                                                                                      (5)

где:  – максимальное значение дисперсии из двух образцов УАБ;

 – минимальное значение дисперсии.

На вход данного блока подаются дисперсии двух сравниваемых образцов УАБ с точностью до 4 знака после запятой, выходом является значение критерия Фишера, что и описывается в блоке «Язык программирования» (здесь также определяется максимальное и минимальное значение дисперсий):

input dx1, dx2;

if dx1 > dx2 then F = dx1/dx2 else F = dx2/dx1;

output F;

В случае равенства дисперсий величина F имеет стандартное распределение Фишера (F-распределение) с числом степеней свободы (n–1, m–1) и доверительной вероятностью P_д = 0.95. Таким образом, критерий F должен лежать в диапазоне значений:

.

Если гипотеза о равенстве дисперсий подтверждена, то на выходе модели будет получено логическое значение 1, что соответствует тому, что сравниваемые образцы УАБ могут быть далее сопоставлены по степени автономности полёта при помощи критерия Стьюдента.

В таблице 5 представлены значения критерия Фишера для каждой пары сравниваемых образцов УАБ.

Таблица 5.

Значение критерия Фишера при парном сравнении образцов УАБ

Анализируя данную таблицу, видно, что значения критерия Фишера не для всех пар сравниваемых УАБ принадлежат допустимому диапазону значений.

Дальнейшее сопоставление образцов УАБ по степени автономности с использованием критерия Стьюдента будет проводиться только для пар, удовлетворяющих условию попадания значений критерия Фишера в допустимый диапазон.

Для этого была построена схема автоматики «Сопоставление образцов УАБ по степени автономности полёта с использованием критерия Стьюдента», которая показана на рисунке 6:

Рисунок 6. Схема «Сопоставление образцов УАБ по степени автономности полёта с использованием критерия Стьюдента»

Скрипт проект, соответствующий данной схеме, содержит код:

h = 1e-3; // шаг интегрирования

//

X1 = [0.63, 0.47, 0.39, 0.44, 0.74, 0.36, 0.40, 0.49, 0.41, 0.41, 0.52, 0.67, 0.45, 0.34, 0.34, 0.50, 0.49, 0.43, 0.44, 0.55, 0.27, 0.59, 0.43, 0.54, 0.42, 0.36, 0.52, 0.60, 0.51, 0.52, 0.54, 0.35, 0.71, 0.77, 0.69, 0.52, 0.59, 0.42, 0.61, 0.51, 0.31, 0.46, 0.51, 0.49, 0.52, 0.58, 0.45, 0.57, 0.51, 0.65];

X3 = [0.29, 0.33, 0.29, 0.26, 0.27, 0.31, 0.30, 0.32, 0.25, 0.24, 0.27, 0.22, 0.26, 0.35, 0.35, 0.33, 0.27, 0.22, 0.32, 0.32, 0.30, 0.33, 0.30, 0.31, 0.28, 0.29, 0.40, 0.30, 0.39, 0.33, 0.32, 0.26, 0.30, 0.25, 0.32, 0.30, 0.32, 0.26, 0.34, 0.31, 0.30, 0.29, 0.36, 0.32, 0.17, 0.26, 0.33, 0.29, 0.22, 0.32];

X8 = [0.54, 0.53, 0.51, 0.56, 0.49, 0.48, 0.48, 0.55, 0.47, 0.50, 0.51, 0.52, 0.57, 0.49, 0.41, 0.54, 0.50, 0.49, 0.51, 0.47, 0.52, 0.47, 0.49, 0.49, 0.53, 0.51, 0.54, 0.50, 0.39, 0.55, 0.47, 0.47, 0.51, 0.57, 0.52, 0.46, 0.48, 0.50, 0.50, 0.41, 0.50, 0.50, 0.48, 0.47, 0.42, 0.42, 0.46, 0.50, 0.53, 0.45];

X9 = [0.19, 0.41, 0.47, 0.38, 0.27, 0.46, 0.22, 0.33, 0.47, 0.37, 0.35, 0.27, 0.65, 0.12, 0.33, 0.52, 0.50, 0.44, 0.40, 0.33, 0.50, 0.36, 0.47, 0.56, 0.38, 0.54, 0.42, 0.44, 0.37, 0.42, 0.33, 0.56, 0.30, 0.16, 0.31, 0.16, 0.44, 0.42, 0.45, 0.37, 0.41, 0.44, 0.49, 0.50, 0.51, 0.29, 0.39, 0.24, 0.44, 0.40];

//

// Pd - доверительная вероятность = 0.95

tmin = -1.9845; // квантиль критерия Стьюдента ((1-Pd)/2, n+m-2)

tmax = 1.9845; // квантиль критерия Стьюдента ((1+Pd)/2, n+m-2)

Данная схема представляет собой реализацию формулы критерия Стьюдента [3, с. 73]:

                                                                                                               (6)

где: ,  – математического ожидания и дисперсия первого сравниваемого образца УАБ;

,  –математического ожидания и дисперсия второго сравниваемого образца УАБ;

n, m – количество элементов в первой и второй выборке.

Критерий Стьюдента, как и критерий Фишера, имеет ограничения на диапазон допустимых значений. Будем считать сопоставление по автономности полёта двух образцов УАБ выполненным, если критерий Фишера, принятый при значении доверительной вероятности P_д = 0.95, с числом степеней свободы (n+m–2) принадлежит интервалу: .

Для отобранных пар образцов 1 и 9, 3 и 8 были рассчитаны значения критерия Стьюдента и сделан вывод о сопоставимости по автономности полёта образцов УАБ. Результаты представлены в таблице 6.

Таблица 6.

Значение критерия Стьюдента при парном сравнении образцов УАБ

В результате видим, что ни одна из сравниваемых пар образцов УАБ не удовлетворяет указанным условиям, потому что значения критерия Стьюдента не попадают в допустимый диапазон значений.

Можно сделать вывод, что сравниваемые образцы УАБ имеют различия по степени автономности полёта и их технические характеристики нельзя считать сопоставимыми. Результаты, отраженные в данной работе, не отличаются от тех, что были бы получены в одном из вышеупомянутых статистических пакетов – Statistica, что свидетельствует о том, что среда SimInTech удобна для решения задач обработки экспериментальных данных и даёт отличную возможность пользователям при должном уровне теоретических знаний погрузиться в процесс решения задачи.

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

1.  Карташов Б.А. Среда динамического моделирования технических систем SimInTech: Практикум по моделированию систем автоматического регулирования: учебное пособие / Б.А. Карташов, Е.А. Шабаев, О.С. Козлов, А.М. Щекатуров – М.: ДМК Пресс, 2017. – 424 с.
2. Справочная система SimInTech [Электронный ресурс]. // SimInTech. URL: http://help.simintech.ru/ (дата обращения: 10.12.2019).
3. Евдокименков В.Н., Базлев Д.А. Автоматизация обработки результатов экспериментальных исследований бортовых комплексов ЛА: Учебное пособие. М.: Изд. МАИ-ПРИНТ, 2011. – 229 с.: ил.