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

Ax +x +Aх b +.

Вычитая из полученного неравенства уравнение для точной системы Ax = b (1.26), получим:

x +Aх .

Перенося слагаемые и умножая обе части на обратную матрицу ( A)^-1 , получаем искомое выражение для х:

х = ( A)^-1( - x ).(2.3)

Для оценки общей величины вектора погрешностейх предложено рассматривать его евклидову норму:

х= (х₁²+ х₂²+…+х_n²)^0,5.

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

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

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

Отсюда вытекает простейший способ уменьшениях путем непосредственного уменьшения погрешностей и .

2.2 Линейное симметричное сглаживание исходных данных расчетной модели. Формулы для неравномерной сетки

Допустим, дискретные медленно изменяющиеся опытные данные {f(t_i)} = {f(t₁), f(t₂), …, f(t_i),…} содержат высокочастотные шумы с неизвестными статистическими параметрами. Получаемые опытные данные обычно являются результатом измерений и предпроцессорной обработки характеристик некоторого реального непрерывного процесса.

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

В классическом варианте сетка значений времени t₁, t₂, …, t_i,… является равномерной (при рассмотренном выше управлении с прогнозированием сетка существенно неравномерная, поскольку регистрируются моменты выполнения некоторых событий).

Так как новые значения сигнала {g(t_i)} рассчитываются непосредственно по прежним значениям {f(t_i)} без перехода в частотную область, то сглаживание относится к пространственным (временным) методам фильтрации.

Самым простым методом сглаживания является скользящее среднее (Moving average filtering), в котором каждый член ряда заменяется простым или взвешенным средним n соседних членов, где n - ширина "окна".

Простое среднее (Simple moving average, SMA - простое скользящее среднее ) является арифметическим средним, определяемым по нечетному числу n сигналов, симметрично расположенных относительно текущей точки. Обозначив k = [n/2], формулу для сглаживания по простому среднему можно представить в виде:

g(t_i) = [f(t_i-k)+…+ f(t_i)+… f(t_i+k)]/n. (2.4)

При простом сглаживании все значения сигнала f(t_i-k),…,f(t_i+k), участвующие в расчете g(t_i), имеют одинаковые коэффициенты, равные 1/ n.

Положительным качеством данного метода является простота и минимальное число вычислений, поскольку при расчетах можно использовать рекурсивное соотношение: g(t_i) = g(t_i-1) + [f(t_i+k) - f(t_i-k-1)]/n.

Главным недостатком простого среднего является то, что на новое значение сигнала g(t_i) в равной степени влияют все прежние сигналы, как близкие к моменту t_i, так и удаленные от него.

Для того, чтобы увеличить вклад близких значений и уменьшить - для удаленных, в формуле для g(t_i) используют весовые коэффициенты, значения которые задают в зависимости от удаленности от средней точки по линейному закону (Linear Weighted Moving Average, LWMA - линейно-взвешенное скользящее среднее), экспоненте (exponential moving average, EMA - экспоненциальное скользящее среднее ), синусоиде и т.д.

Рассмотрим общий случай линейно-взвешенного закона усреднения. Если принять относительный коэффициент при центральном i-том сигнале через 1, а относительную величину убывания коэффициента при сдвиге на один номер узла обозначить через h, то с учетом n=2k+1 сумма всех коэффициентов при узлах равна:

S = 1+2[(1-h)+(1-2h)+…+(1-kh)] = (1+2k) - hk(k+1) = n-hk(k+1).

При этом формула для сглаженного значения сигнала принимает вид:

g(t_i)=[f(t_i)+(1-h)(f(t_i-1)+(t_i+1))+(1-2h)(f(t_i-2)+(t_i+2))+…+(1-kh)(f(t_i-k)+ f(t_i-+k)]/S. (2.5)

В частном случае при 1- kh = h; h =1/(k+1) формула принимает вид:

S = n - k = k+1; g(t_i) =[f(t_i)+ k(f(t_i-1)+(t_i+1))+(k-1)(f(t_i-2)+(t_i+2))+ …+ 1(f(t_i-k)+ f(t_i-+k)))]/( k+1).(2.6)

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

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

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

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

Рассмотрим симметричную относительно центральной точки t_j линейно-взвешенную фильтрацию по общему нечетному числу m=2s+1 сигналов.

Обозначим через _{i j} отклонение текущей точки t_i относительно центральной точки t_j :

_{i j} = t_i - t_j,(2.7)

а через _m - максимальное отклонение текущих точек t_i относительно центральной точки t_j .

Очевидно, для _m есть две возможности:

1) _m = _{(j+s) j} - максимум достигается в точке t_j = t_j+s (рис.2.2),

2) _m = _{(j-s) j} - максимум достигается в точке t_j = t_j-s (рис.2.3).

Рис.2.2. Случай 1) _m = _{(j+s) j}

Рис.2.3. Случай 2) _m = _{(j-s) j}

Как в первом, так и во втором случае линейный коэффициент для узлов рассчитывается по общей формуле:

k_i = 1 - _{i j} / _m.(2.8)

Рассмотрим случай 1). Учитывая равенство линейного коэффициента 1 в центральном узле и 0 - в самом отдаленном, формулу для сглаженного значения сигнала можно выразить через линейные коэффициенты следующим образом:

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

(2.9)

Аналогично выводим формулу для случая 2):

(2.10)

Полученные формулы (2.9) и (2.10) требуют выполнения минимального числа элементарных расчетных операций. Суммируя числа операций, получаем, что при выполнении симметричной относительно центральной точки линейно-взвешенной фильтрации минимальные числа операций в одной расчетной точек следующие:

- сложение и вычитание: s(s) = 4s-2;

- умножений: m(s)=2s;

- делений: d(s)=1.

2.3 Основные допущения при построении имитационной модели

Симметричная линейно-взвешенная фильтрация, описанная выше, имеет единственный параметр s, задающий количество информативных точек t_i, лежащих по каждую сторону от центрального узла t_j. Очевидно, уменьшение погрешности в препроцессорных данных и существенно зависит от величины s.

В доступных источниках информации нет данных о том, каким образом уменьшается погрешность входных данных в зависимости от:

1) полуширины фильтра s,

2) типа базовой функции f(t), на которую накладывается погрешность,

2) равномерноcть или неравномерность сетки значений узлов по времени t₁, t₂, …, t_i, среднее значение ее шага, при неравномерной сетке - дополнительные статистические ее характеристики,

3) параметров самой погрешности f(t) в узлах t₁, t₂, …, t_i.

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

I. Рассмотрены значения полуширины фильтра s, равные s =1,2,3.

II. Для характеристики изменения времени взят сегмент [0,10]. На нем в качестве базовых (опорных) функций f(t) приняты:

1) константа f₁(t)0,

2) линейная функция с малой скоростью роста f₂(t)0,1t,

3) линейная функция со средней скоростью роста f₃(t)t,

4) линейная функция с высокой скоростью роста f₄(t)10t,

5) квадратичная функция с малым коэффициентом роста f₅(t) 0,001t²,

6) квадратичная функция со средним коэффициентом роста f₆(t) 0,01t²,

7) квадратичная функция со значительным коэффициентом роста f₇(t) 0,1t²,

8) квадратичная функция с большим коэффициентом роста f₈(t) t²,

9) синусоидальная функция с малой осцилляцией f₉(t)sin(0,01t),

10) синусоидальная функция со средней осцилляцией f₁₀(t) sin(0,1t),

11) синусоидальная функция с высокой осцилляцией f₁₁(t)sint.

III, Для того, чтобы учесть при выборе метода моделирования моментов наступления событий в реальной системе инертность функционирования самой системы, ограничения на минимальную и максимальные скорости выполнения процессов в ней, предложены следующие соотношения.

1. Масштабное задание начального нулевого значения времени события t₀ = 0 и временного отрезка t₁ = t₁ - t₀ = 1 между нулевым и первым событиями.

2. Соотношение между последующими отрезками между событиями j и предыдущим j-1: t_j = (0,7+0,6)t _{j -1}, где - случайная величина, равномерно распределенная на отрезке [0,1], j >1 .

3. Ограничения на минимальную и максимальную скорости процессов, происходящих с объектом управления.

Если t_j < 0,33, то принимаем: t_j =0,33.

Если t_j >3, то принимаем: t_j =3.

Погрешности исходных данных - нормально распределенные случайные величины. Их значения рассчитываются по второму варианту метода Бокса-Мюллера.

2.4 Алгоритм имитационного моделирования зашумленных данных и их симметричного линейного сглаживания. Анализ результатов

1. Всего рассмотрено 11 опорных функций, перечисленных в п.2.3, для которых принят масштабный фактор, равный 1.

2. По каждой функции для каждого значения масштабного коэффициента шумов СОЕ (0,1 - высокие шумы, 0.01 - средние шумы, 0,001 - малые шумы) задано 10 испытаний.

3. В каждом испытании моделируется опорная функция и идеальная шумовая составляющая:

а) вначале строится 50 значений моментов времени t[i](i=0,…,49),

б) для каждого момента t[i] рассчитывается точное значение опорной функции f_k(t[i]), шумовая составляющая f_k(t[i]) и полная величина сигнала fv_k(t[i]) = f_k(t[i]) + f_k(t[i]).

4. В каждом испытании:

а) для 30 значений времени t[i](i=11,…,40) моделируется процесс линейной симметричной фильтрации со значениями полуширины фильтра s=1,2,3, при этом получаются массивы значений f1(t[i]), f2(t[i]) , f3(t[i]) (i=11,…,40),

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

S_E,S_E1, S_E2, S_E3, а также средние значения среднеквадратичного отклонения исходной зашумленной функции S_SG и относительные доли уменьшения данной величины после сглаживания, которые рассчитывались по формулам:

DS1=(S_SG-S_SG1)/S_SG;

DS2=(S_SG-S_SG2)/S_SG;

DS3=(S_SG-S_SG3)/S_SG.

Отрицательные значения DS1,DS2,DS3 свидетельствуют о том, что после фильтрации среднеквадратическое отклонение получаемой кривой от опорной функции увеличилось, т.е. качество данных обратной связи в результате фильтрации по заданному методу ухудшилось.

Сводные результаты моделирования даны в таблице 2.1.

Таблица 2.1

№ фнк	СОE	S_E	S_E1	S_E2	S_E3	S_SG	DS1	DS2	DS3=
1	0,1	0.289	0.297	0.301	0.293	0.317	0.048	0.300	0.480
	0,01	0.029	0.030	0.031	0.030	0.032	0.048	0.300	0.480
	0,001	0.003	0.003	0.004	0.003	0.003	0.048	0.300	0.480
2	0.1	0.473	0.481	0.470	0.482	1.061	0.096	0.489	0.603
	0.01	0.047	0.048	0.047	0.047	0,106	0.040	0.446	0.536
	0,001	0.005	0.005	0.004	0.005	0.0106	-1.24	-1.23	-1.65
3	0,1	0.257	0.271	0.277	0.289	2.578	0.057	0.359	0.420
	0,01	0.026	0.040	0.051	0.063	0.026	-1.62	-1.94	-2.21
	0,001	0.003	0.016	0.028	0.041	0.003	-22.9	-26.8	-30.02
4	0,1	0.420	0.584	0.630	0.829	0.564	-2.45	-2.47	-3.16
	0,01	0.042	0.219	0.278	0.481	0.056	-30.1	-31.7	-39.7
	0,001	0.004	0.182	0.243	0.446	0.006	-307	-324	-405
5	0,1	0.195	0.194	0.200	0.204	0.171	0.102	0.386	0.519
	0,01	0.020	0.019	0.020	0.022	0.017	0.069	0.358	0.494
	0,001	0.002	0.002	0.002	0.003	0.002	-3.16	-2.30	-2.78
6	0,1	0.245	0.250	0.257	0.261	0.303	0.086	0.500	0.651
	0,01	0.024	0.027	0.030	0.034	0.030	0.123	0.560	0.535
	0,001	0.002	0.005	0.007	0.011	0.003	-2.94	-3.91	-5.44
7	0,1	0.277	0.260	0.264	0.292	0,375	-0.07	0.18	0.175
	0,01	0.028	0. 01	.0298	0.061	.0375	-3.60	-4.79	-6.09
	0,001	0.003	-0.01	0.007	0.038	0.004	-43.1	-56.3	-69.7
8	0,1	0.307	0.672	1.089	1.645	.496	-3.98	-4.35	-5.70
	0,01	0.031	0.385	0.81	1.36	0.050	-46.5	-50.6	-64.2
	0,001	0.003	0.356	0.781	1.33	0.005	-472	-512	-649
9	0,1	0.190	0.201	0.211	0.214	0.230	0.054	0.429	0.623
	0,01	0.019	0.020	0.02	0.021	0.023	0.057	0.427	0.623
	0,001	0.002	0.002	0.002	0.002	0.002	-0.02	0.208	0.207
10	0,1	0.331	0.337	0.333	0.333	0.451	0.036	0.417	0.552
	0,01	0.033	0.034	0.033	0.032	0.045	0.037	0.421	0.557
	0,001	0.003	0.003	0.003	0.002	0.005	0.040	0.428	0.497
11	0,1	0.360	0.362	0.367	3.64	0.464	0.024	0.371	0.478
	0,01	0.036	0.035	0.034	0.030	0.047	-0.05	-0.15	-0.97
	0,001	0.004	0.003	0.0002	-0.003	0.005	-2.87	-7.73	-17.3

Анализ результатов выполненного моделирования показывает следующее.

1. Константа f₁(t)0. При увеличении s средняя величина шума Е мала и практически не изменяется. Среднее отклонение монотонно уменьшается на 5%, 30%, 48%.

2. Линейная функция с малой скоростью роста f₂(t)0,1t . Картина аналогична константе при больших и средних погрешностях. Для малых значений среднее отклонение ухудшается примерно в 2 раза.

3. Линейная функция со средней скоростью роста f₃(t)t. Картина подобна константе только при большой погрешности . При средней с ростом s показатели ухудшаются незначительно. При малой показатели значительно ухудшаются, особенно по - она увеличивается в 22-30 раз.

4. Линейная функция с высокой скоростью роста f₄(t)10t.Результаты неудовлетворительны при всех величинах погрешностей и s=1,2,3.

5. Квадратичная функция с малым коэффициентом роста f₅(t) 0,001t². Результаты удовлетворительны только при больших и средних . При малых в несколько раз вырастает .

6. Квадратичная функция со средним коэффициентом роста f₆(t) 0,01t². Результаты аналогичны предыдущему типу.

7.Квадратичная функция со значительным коэффициентом роста f₇(t) 0,1t². Для большой погрешности сглаживание практически не меняет качественных показателей шумов. При средних и высоких - ухудшает их.

8. Квадратичная функция с большим коэффициентом роста f₈(t) t². Сглаживание ухудшает шумы при всех величинах погрешности.

9. Синусоидальная функция с малой осцилляцией f₉(t)sin(0,01t). Картина аналогична константе.

10. Синусоидальная функция со средней осцилляцией f₁₀(t) sin(0,1t). Картина аналогична константе.

11. Синусоидальная функция с высокой осцилляцией f₁₁(t)sint. При больших погрешностях картина аналогична константе, при средних и малых показатели ухудшаются.

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

Предложено формульное задание для симметричной относительно центральной точки t_j линейно-взвешенной фильтрации по общему нечетному числу m=2s+1 сигналов для случая неравномерной сетки по времени, учитывающее два возможных случая отклонение текущих точек относительно центральной точки. Оценена вычислительная сложность фильтров. Для фильтрации одной точки при полуширине фильтра s требуется: сложений и вычитаний: s(s) = 4s-2; умножений: m(s)=2s; делений: d(s)=1. Т.о. предложенный алгоритм фильтрации имеет линейную сложность по s.

Разработан алгоритм имитационного моделирования зашумленных данных и их симметричного линейного сглаживания. По выполненному алгоритму разработана расчетная программа на языке С++ и выполнены расчеты для 11 видов функций.

Анализ результатов моделирования в целом показал, что:

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

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

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

3. Аппроксимационное сглаживание данных обратной связи

3.1 Постановка задачи аппроксимационного сглаживания

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

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

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

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

Рассмотрим текущий момент времени t_j и построим квадратичную параболу, которая точно проходит через уже сглаженные ранее предыдущие точки траектории, соответствующие моментам времени t_j-1, t_j-2, t_j-3. Число предыдущих точек принято равным 3 по той причине, что у квадратичной параболы 3 коэффициента. В каноническом виде ее уравнение представляют с помощью постоянных коэффициентов в виде следующей формулы:

Р₂(t) = C₀+C₁t+C₂t²,(3.1)

где {C₀, C₁, C₂} - постоянные коэффициенты, C₂ - квадратичный, C₁- линейный, C₀- свободный.

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

Р₂(t) = (C₂t + C₁)t + C₀ (3.2)

на 1 расчет значения квадратичной параболы требуется выполнение:

- сложение и вычитание: s = 2;

- умножений: m=2;

- делений: d=0.

3.2 Вывод основных соотношений квадратичной аппроксимации

Условия точного прохождения параболы через заданные точки с координатами (t_j-1, f_j-1),(t_j-2, f_j-2), (t_j-3, f_j-3) имеют вид:

f _j-1 = f (t_j-1) = C₀+C₁ t_j-1+C₂ t_j-1²,

f _j-2 = (t_j-2) = C₀+C₁ t_j-2+C₂ t_j-2²,

f _j-3 = f (t_j-3) = C₀+C₁ t_j-3+C₂ t_j-3². (3.3)

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

f _j-1 = C₁t_j-1+C₂ t_j-1(t_j-1+ t_j-2),

f _j-2 = C₁t_j-2+C₂ t_j-2(t_j-2+ t_j-3), (3.4)

где t_j-1= t_j-1 - t_j-2 , t_j-2= t_j-2 - t_j-3 , f _j-1= f _j-1 - f _j-2 , f _j-2= f _j-2 - f _j-3 .

Умножая первое уравнение на t_j-2, а второе - на t_j-1, приводим систему (3.4) к виду:

f _j-1t_j-2 = C₁t_j-1t_j-2+C₂ t_j-1t_j-2 (t_j-1+ t_j-2),

f _j-2t_j-1 = C₁t_j-1t_j-2+C₂ t_j-1t_j-2(t_j-2+ t_j-3). (3.5)

Вычитая из второго уравнения первое, получим:

f _j-1t_j-2 - f _j-2t_j-1 = C₂ t_j-1t_j-2 (t_j-3- t_j-1).

Поскольку все сомножители при C₂ не равны 0, то решение полученного уравнения относительно C₂ всегда существует и равно:

C₂ =(f _j-1t_j-2 - f _j-2t_j-1)/ t_j-1t_j-2 (t_j-3- t_j-1).(3.6)

Величину линейного коэффициента C₁ по уже известному значению C₂ с наименьшими вычислительными затратами можно определить из первого уравнения системы (3.4):

C₁ = f _j-1/t_j-1 - C₂(t_j-1+ t_j-2) .(3.7)

По найденным коэффициентам C₂ и C₁ свободный коэффициент C₀ находим из первого уравнения системы (3.3):

C₀ = f _j-1 - C₁ t_j-1 - C₂ t_j-1². (3.8)

Переходя в расчете (3.8) к формуле Горнера, получим расчетную схему вычисления коэффициентов аппроксимирующей параболы с наименьшим числом расчетных операций:

t_j-1= t_j-1 - t_j-2 , t_j-2= t_j-2 - t_j-3 , f _j-1= f _j-1 - f _j-2 , f _j-2= f _j-2 - f _j-3,

C₂ =(f _j-1t_j-2 - f _j-2t_j-1)/ t_j-1t_j-2 (t_j-3- t_j-1),

C₁ = f _j-1/t_j-1 - C₂(t_j-1+ t_j-2),

C₀ = f _j-1 - (C₂ t_j-1 + C₁) t_j-1.(3.9)

Полученные соотношения (3.9) позволяют рассчитать коэффициенты аппроксимирующей параболы {C₀, C₁, C₂}по значениям уже сглаженных значений траектории, соответствующих моментам времени t_j-1, t_j-2, t_j-3.

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

- сложение и вычитание: s = 10;

- умножений: m=7;

- делений: d=2.

3.3 Варианты алгоритмов квадратичного аппроксимационного сглаживания данных

Рассмотрим текущее значение времени t_j, в котором соответствующее значение функции f_j содержит шумы. Необходимо найти его значение g_j, которое в максимально возможной степени очищено от высокочастотных шумов.

Для всех вариантов аппроксимационных квадратичных алгоритмов вначале по предыдущим моментам времени t_j-1, t_j-2, t_j-3 рассчитываются коэффициенты {C₀, C₁, C₂} аппроксимирующей параболы Р₂(t) и определяется значение данной параболы Р₂(t_j) в текущей точке t = t_j.

В зависимости от использования в фильтрации последующих точек {t_j+1, t_j+2, t_j+3 }, в работе предложено 4 варианта метода, эффективность которых исследована на тех же модельных кривых, что и в случае линейной фильтрации. Число используемых при сглаживании последующих точек назовем глубиной метода и обозначим ее h.

Также свойства алгоритма зависят от того, с какой долей в общую сумму в слагаемых центральной и последующих точках входят значения самой функции f_i и ее значение Р₂(t_i), аппроксимированное по параболе . Обозначим долю аппроксимированного значения Р₂(t_i) через , а долю исходного значения f_i примем равной (2-) для того, чтобы в сумме они давали 2 и величина нормирующего знаменателя формулы не зависела от .

Теоретически величина может изменяться на отрезке значений от 0 до 2. При =0 аппроксимация не учитывается и сглаживание становится асимметричным по последующим значениям. При =2 полностью исключаются исходные значения функции и получаемые значения теряют смысл. Определение оптимальной доли аппроксимированного значения Р₂(t_i) исследовано в следующем параграфе.

Рассмотрим семейство аппроксимационных квадратичных методов с глубинами h =0,1,2,3 и заданной долей аппроксимированного значения .

I. h = 0. При нулевой глубине последующие точки не учитываются при фильтрации в центральной точке.

1. После выполнения общих предварительных операций искомое значение g_j рассчитывается по формуле:

g_j⁰ = ((2-)f_j + Р₂(t_j))/2. (3.10)

Для учета последующих точек при сглаживании в центральной точке t=t_j вначале по предыдущим точкам рассчитывается полуширина полосы фильтра по точкам аппроксимации: t_m = t_j - t_j-3 .

II. h = 1. При сглаживании учитывается только одна последующая точка t_{j +1}, соседняя с центральной. После выполнения предварительных действий производятся следующие расчеты.

1. Расчет линейного t_{j +1}весового коэффициента при первой последующей точке t = t_{j +1}.: k₁ = (1-(t_{j +1} - t_j)/ t_m).

2. Расчет значения параболы Р₂(t_j+1) в точке t = t_{j +1}.

3. Расчет значения g_j по формуле:

g_j ¹ = 0,5((2-)f_j + Р₂(t_j) + k₁ ((2-)f_j+1+ Р₂(t_j+1)))/(1+ k₁). (3.11)

III. h = 2. При сглаживании учитываются две последующие точки. Строится аналогично метолу с одной точкой.

1.Расчет весовых коэффициентов для двух последующих точек t = t_j+1 и t = t_{j +2}.: k₁ = (1-(t_{j +1} - t_j)/ t_m), k₂ = (1-(t_{j +2} - t_j)/ t_m).

2. Расчет значений параболы Р₂(t_j+1), Р₂(t_j+2) в точках t = t_{j +1}, t = t_{j +2}.

3. Расчет значения g_j по формуле:

g_j² = 0,5((2-)f_j+Р₂(t_j) + k₁((2-)f_j+1+Р₂(t_j+1)) + k₂ ((2-)f_j+2+ Р₂(t_j+2)))/ (1+k₁+k₂). (3.12)

IV. h = 3. Аппроксимационный квадратичный метод с использованием трех последующих точек.

1.Расчет весовых коэффициентов для трех последующих точек от t = t_j+1 до t = t_{j +3}.: k₁ = (1-(t_{j +1} - t_j)/ t_m), k₂ = (1-(t_{j +2} - t_j)/ t_m), k₃ = (1-(t_{j +3} - t_j)/ t_m).

2. Расчет значений параболы Р₂(t_j+1), Р₂(t_j+2), Р₂(t_j+3), в точках t = t_{j +1}, t = t_{j +2}, t = t_{j +2}.

3. Расчет значения g_j по формуле:

g_j³ = 0,5((2-)f_j+Р₂(t_j) + k₁((2-)f_j+1+Р₂(t_j+1)) + k₂ ((2-)f_j+2+Р₂(t_j+2)) + k₃ ((2-)f_j+3+Р₂(t_j+3)))/(1+k₁+k₂+k₃). (3.13)

При сравнительном расчете всех сглаженных значений g_j⁰, g_j¹, g_j², g_j³ в расчетной программе был использован сокращенный алгоритм расчета, исключающий повторные включения в расчеты линейных весовых коэффициентов и значений функции.

Обозначая через nc0 - nс3 нормирующие коэффициенты рассматриваемых 4 методов, а через dc - линейный коэффициент текущей присоединяемой последующей точки, получаем следующие сокращенные расчетные схемы.

0. h = 0. dc=1; nco=1;

g_j⁰=0,5((2-)f_j + Р₂(t_j)).

1. h = 1. dc=1-( t_{j +1} - t_j )/dtm; nc1=nc0+dc;

g_j¹=( g_j⁰ +0.5*dc*((2-)f_j+2+ Р₂(t_j+2)))/ nc1.

2. h = 2. dc=1-( t_{j +2} - t_j )/dtm; nc2=nc1+dc;

g_j²=( g_j¹ +0.5*dc*((2-)f_j+2+ Р₂(t_j+2)))/ nc2.

3. h = 3. dc=1- ( t_{j +3} - t_j )/dtm; nc3=nc2+dc;

g_j³=( g_j²+0.5*dc*((2-)f_j+3+ Р₂(t_j+3)))/ nc3 (3.14)

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

1. h = 0. При выполнении предварительных расчетов и итоговом расчете (3.10) общее число операций, затрачиваемых на одно квадратичное сглаживание при нулевой глубине следующее:

- сложение и вычитание: s = 14;

- умножений: m=12;

- делений: d=3.

2. h = 1. При глубине, равной 1, общее число элементарных операций, затрачиваемых на одно квадратичное сглаживание:

- сложение и вычитание: s = 20;

- умножений: m=14;

- делений: d=5.

3. h = 2:

- сложение и вычитание: s = 25;

- умножений: m=16;

- делений: d=6.

4. h = 3

- сложение и вычитание: s = 32;

- умножений: m=17;

- делений: d=7.

3.4 Алгоритм имитационного моделирования зашумленных данных при помощи квадратичного аппроксимационного сглаживания. Анализ результатов

Исследуем на примере постоянной функции f(t)=0 влияние доли аппроксимированного значения на качество сглаживания. Для оценки качества сглаживания по результатам расчета математических ожиданий шумов Е, Е0, Е1, Е2, Е3 (для исходных нормальных шумов и остаточных шумов после сглаживания значений функции на глубину h=0,1,2,3) были выбраны относительное величины уменьшения математических ожиданий шумов после сглаживания:

DЕ0=(Е0-Е)/Е;

DЕ1=(Е1-Е)/Е;

DЕ2=(Е2-Е)/Е;

DЕ3=(Е3-Е)/Е .

В качестве интегральных критериев качества приняты:

1) суммарное уменьшение матожиданий шумов при фильтрации всей совокупностью методов (h=0,1,2,3)

SUMD= DЕ0 +DЕ1 + DЕ2 + DЕ3,

2) минимальное значение уменьшения матожиданий шумов при фильтрации всей совокупностью методов (h=0,1,2,3)

MD=min(DЕ0,DЕ1, DЕ2, DЕ3).

Выборочные расчеты при значениях от 0.2 до 1.5 показали, что более информативным является второй критерий. Усредненные значения данного критерия для функции f(t)=0 приведены в таблице 3.1.

Таблица 3.1.Усредненные значения критерия MD для функции f(t)=0

	0,2	0,3	0,4	0,5	0,6	0,7	0,8	0,9	1,0	1,1	1,2	1,3	1,4
MD	0,04	0,03	0,025	0,02	0,015	0,03	0,04	0,04	0,05	0,05	0,06	0,08	0,11

В качестве оптимального было принято значение = 0,6.

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

1. Как и главе 2, рассмотрено 11 опорных функций, перечисленных в п.2.3, для которых принят масштабный фактор, равный 1.

2. По каждой функции для каждого значения масштабного коэффициента шумов СОЕ (0,1 - высокие шумы, 0.01 - средние шумы, 0,001 - малые шумы) задано 10 испытаний.

3. В каждом испытании моделируется опорная функция и идеальная шумовая составляющая:

а) вначале строится 50 значений моментов времени t[i](i=0,…,49),

б) для каждого момента t[i] рассчитывается точное значение опорной функции f_k(t[i]), шумовая составляющая f_k(t[i]) и полная величина сигнала fv_k(t[i]) = f_k(t[i]) + f_k(t[i]).

4. В каждом испытании:

а) для 30 значений времени t[i](i=11,…,40) моделируется процесс квадратичного аппроксимационного сглаживания при значениях глубины h=0,1,2,3, при этом получаются массивы значений f0(t[i]), f1(t[i]), f2(t[i]) , f3(t[i]) (i=11,…,40),

б) по результатам моделирования исходной зашумленной функции и ее сглаженных вариантов определяются усредненные математические ожидания S_E, S_E0, S_E1, S_E2, S_E3, по ним дополнительно рассчитываются относительные коэффициенты снижения среднего уровня шума

DЕ0=(Е0-Е)/Е;

DЕ1=(Е1-Е)/Е;

DЕ2=(Е2-Е)/Е;

DЕ3=(Е3-Е)/Е;

а также средние значения среднеквадратичного отклонения исходной зашумленной функции S_SG и относительные доли уменьшения данной величины после сглаживания, которые рассчитывались по формулам:

DS0=(S_SG-S_SG0)/S_SG;

DS1=(S_SG-S_SG1)/S_SG;

DS2=(S_SG-S_SG2)/S_SG;

DS3=(S_SG-S_SG3)/S_SG.

Получение отрицательных значений у относительных значений DS0,DS1,DS2,DS3 свидетельствуют о том, что после сглаживания среднеквадратическое отклонение получаемой кривой от опорной функции увеличилось, т.е. разброс значений увеличился.

Текст программы на языке С++, в котором по сокращенному алгоритму (3.14) рассчитываются сглаженные значения 11 модельных функций, а также рассчитываются их математические ожидания и доли уменьшения их среднеквадратичных отклонений после их квадратичного аппроксимационного сглаживания с глубинами h=0,1,2,3.

Сокращенные результаты имитационного моделирования аппроксимационных квадратичных алгоритмов сглаживания данных приведены в таблице 3.2.

Таблица 3.2 Сводные результаты моделирования аппроксимационных квадратичных алгоритмов сглаживания данных

№ фнк	СОE	S_E	DS_E0	DS_E1	DS_E2	DS_E3
		S_SG	DS0	DS1	DS2	DS3
1	0,1	0.206	0.206	0.0758	0.1449	0.471
		0.406	-0.600	-1.295	-1.1675	-0.231
	0,01	0.031	0.048	0.545	18.2943	37.05
		0.041	-0.596	-1.292	-2.38797	-4.830
	0,001	0.003	0.0480	5.019	180.39	365.08
		0.004	-0.595	-1.503	-26.03	-56.196
2	0.1	0.233	0.015	0.119	= 3.519	7.4167
		0.271	-0.725	-1.327	-1.372	-1.334
	0.01	0.023	0.015	1.366	32.00	68.42
		0.002	0.027	-0.725	-1.376	-8.309
	0,001	= 0.015	0.015	13.831	316.82	678.49
		0.003	-0.725	-2.989	-87.809	-188.40
3	0,1	0.203	0.009	1.412	37.402	78.608
		0.125	-0.372	-0.662	-15.105	-37.562
	0,01	0.020	0.009	13.762	371.31	779.77
		0.012	-0.372	-3.671	-162.93	-385.38
	0,001	0.002	0.009	137.27	3710.38	7791.4
		0.001	-0.371	-44.69	-1641.7	-3863.7
4	0,1	0.203	0.009	1.412	37.402	78.608
		0.125	-0.372	-0.662	-15.105	-37.562
	0,01	0.020	0.009	13.762	371.31	779.77
		0.001	-0.371	-44.69	-1641.7	-3863.7
	0,001	0.260	0.009	137.279	3710.4	= 7791.4
		0.001	-0.371	-44.689	-1641.7	-3863.7
5	0,1	0.260	0.104	0.534	3.8806	8.800
		0.284	-0.556	-1.335	-2.0689	-3.830
	0,01	0.026	0.104	4.063	36.992	81.78
		0.028	-0.556	-1.698	-22.922	-47.93
	0,001	0.006	0.104	39.352	368.10	811.54
		0.003	-0.556	-12.096	-239.05	-490.1
6	0,1	0.818	0.048	0.223	2.355	4.999
		1.804	-0.365	-0.8450	-0.908	-0.683
	0,01	0818	0.0478	1.676	21.06	44.51
		0.180	-0.365	-0.909	-5.061	-10.13
	0,001	0.008	0.0478	16.208	208.1	439.7
		0.018	-0.365	-2.281	-56.95	-109.6
7	0,1	0.328	0.031	2.312	32.93	68.99
		0.409	-0.254	-0.709	-12.70	-25.67
	0,01	0.033	0.031	22.70	326.78	684.13
		0.041	-0.254	-2.859	-137.1	-266.59
	0,001	0.003	0.031	226.55	3265.2	6835.8
		0.004	-0.254	-34.49	-1382.4	-2675.8
8	0,1	0.258	0.0075	24.87	278.1	618.85
		0.268	-0.528	-7.262	-173.9	-360.04
	0,01	0.026	0.0075	248.9	2778.4	6183.5
		0.027	-0.528	-82.16	-1747.4	-3609.0
	0,001	0.003	0.0075	2489.6	27780.9	61829.6
		0.003	-0.528	-834.14	-17482.3	-36099.0
9	0,1	0.427	0.0519	0.0819	0.44552	1.0178
		0.412	-0.378	-0.900	-0.7452	0.1297
	0,01	0.0427	0.0519	0.167	2.0596	4.5011
		0.042	-0.378	-0.896	-0.8722	-0.6360
	0,001	0.004	0.0520	1.023	18.200	39.334
		0.004	-0.378	-0.878	-6.2412	-13.038
10	0,1	0.305	0.003	0.0681	1.1248	2.0156
		0.262	-0.537	-1.539	-1.334	-0.6956
	0,01	0.031	0.0048	0.513	8.153	15.069
		0.027	-0.538	-1.640	-3.594	-8.0617
	0,001	0.003	0.022	4.967	78.44	145.61
		0.003	-0.542	-5.84	-39.42	-87.72
11	0,1	0.402	0.007	0.0489	0.3660	0.540
		0.612	-0.376	-0.886	-0.8342	-0.273
	0,01	0.040	0.0990	0.652	0.76550	0.0896
		0.061	-0.426	-2.527	-4.80	-7.114
	0,001	0.004	1.158	6.681	4.761	6.387
		0.006	-4.81	-29.9	-55.03	-79.90

Анализ результатов выполненного моделирования показывает следующее.

По всем рассмотренным функциям и масштабам погрешностей картина одинакова: наилучшие результаты дает метод глубины 0. Методы глубины 1, 2, 3 дают либо заметно худшие результаты либо вообще не применимы, поскольку приводят к резкому возрастанию погрешностей.

Уменьшение математического ожидания шумов после применения аппроксимационного квадратического метода глубины 0 составляет от 0.03 у функции №10 (синусоидальная функция со средней осцилляцией f₁₀(t) sin(0,1t)) до 0.104 у функции №5 (линейная функция с высокой скоростью роста f₄(t)10t). Исключение составляет функция №11 (Синусоидальная функция с высокой осцилляцией f₁₁(t)sint) при малых погрешностях, где данный метод привел к небольшому (на 16%) возрастанию математического ожидания шумов.

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

На основе анализа результатов имитационного моделирования дана постановка задачи аппроксимационного сглаживания.