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

УДК 519.254

Энгельгардт В.В.

Самарский государственный университет путей сообщения E-mail: hexware@gmail.com

ЧИСЛЕННЫЙ АНАЛИЗ ГЕНЕТИЧЕСКОГО АЛГОРИТМА В МЕТОДЕ СТРУКТУРНО-ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ ЛИНЕЙНЫХ ДИНАМИЧЕСКИХ СИСТЕМ С ПОМЕХАМИ НА ВХОДЕ И ВЫХОДЕ

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

Ключевые слова: структурно-параметрическая идентификация, линейные динамические системы, эволюционные вычисления, генетический алгоритм, целочисленное программирование.

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

Достаточно сложной задачей в этом случае является определение параметров и структуры модели (порядки по входу и выходу) этого уравнения.

Задача параметрической идентификации динамических систем при наличии помех во входных и выходных сигналах, является более сложной, чем задача регрессионного анализа. В [1] предложен метод для решения таких задач, а рекуррентные модификации разрабатывались и исследовались в [2], [3].

В случае если априорная информация об объекте исследования отсутствует (или она нуждается в подтверждении), существуют методы определения порядка математических моделей динамических систем при наличии помех в выходном сигнале. Например: использование в качестве критерия квадратичной функции потерь [4], [5], [6], [7], информационного критерия Ака-ике (А1С) [8], критерия наименьшего размера модели Риссанена (ЫБЬ)[8], использование матрицы моментов [4], [5], [6], [7]. Можно также использовать информацию о передаточных функциях [6], импульсной характеристике (ИХ) [9, с. 62]. Но данные методы не способны ответить на все вопросы об особенностях системы и обозначают лишь некоторые ориентиры в выборе возможной модели.

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

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

Постановка задачи

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

z-Ib<m>z,-m= Ца^ , (1)

m=1 j=1 m=0

входная и выходная переменная наблюдается с аддитивными помехами в виде:

y. = z () w(j ' = x( > + ^ >(0.

Требуется определить по y, w(J) вектора параметров a, b и порядки г, г..

В качестве функции потерь для оценки параметров ЛРУ в условиях помех использовался следующий критерий [10, c.44]:

min ю-1 (b,a (1)...а (d)bN (b,a (l)...a (а)),где (2) (b}

UN (b,a «...a (d))

Y - AY

b

a

v у

■ Y - Ay

b

a

v >)

(■,) - скалярное произведение;

ra(b,a(1) ...a(d)) =

= Ô12 [l + bTb + y (1) (a(1) )Ta(1) ... y(d) (a(d) )Ta(d) ], ô12 - средняя дисперсия помехи наблюдения ^(г) ;

(ô0))2 - средняя дисперсия помехи наблю-

j ))2

дения )(0, y(j) = ^ ;

О1

a = (a1 ...ad)T , a() = (...aj , b = (b1 ...br)T ,

Y = (1 ...yN)T ,

Уо L yi-r W((1)L W» W((d)L wjdd

Ay,w = MM M : ... M :

У N-1 LYN-r WN (1)LW <N-1 w Nd-1 LW N-rd ■

Приведённый критерий позволяет получить сильно состоятельные оценки [10, с. 43] при выполнении следующих условий:

1. Множество в , которому априорно принадлежат истинные значения параметров устойчивой линейной системы, являются компактом.

2. Помехи ) )(г) j = 1, d статистически независимы и удовлетворяют следующим условиям:

е( (i+1)/ ^ (i0 (i))=0 п.н, где E - оператор математического ожидания; е( (i+1)/ ^ (i0 (i))< h, где h - случайная величина; е(ь)<ж~;

)4(i))<n п.н; е(( +1)/Ç()(ioU(()(i)]= 0 п.н; е(((())2 (i +1)/ Ç(()(i0 )..., Ç(()(i))< h((), h® - случайная величина; E(h(())< п() < E((()(i))< n()1 < ~.

3. {x (1),.x (d)} статистически не зависят от

{()}{{(j\()} j = M .

4. Вектор входных переменных и истинные значения параметров удовлетворяют условиям:

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

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

Задача структурной идентификации формализована таким образом, что её решение закодировано в виде вектора фиксированной длины, где каждый элемент соответствуют сдвигу по r, г.. Таким образом, задача сводится к решению задачи целочисленного программирования (Integerprogramming) которая относятся к классу NP-hard (nondeterministic polynomial time) [12].

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

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

Из полученного множества решений («поколения») с учётом значения «приспособленности» выбираются решения, к которым применяются «скрещивание» и «мутация».

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

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

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

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

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

В общем виде алгоритм структурно-параметрической идентификации имеет следующий вид:

1. Создать базы данных для структур моделей: store - для всех параметрически идентифицированных структур, stab_store - только для устойчивых систем, store_p - является результирующим тезаурусом для моделей с минимальным значением критерия для каждой гра-

d

ницы сложности p ( p = r + ^ r;. );

j=i

2. Задаём ограничения для алгоритма: степень адекватности А, максимально допустимая сложность p;

3. Задаём начальную границу сложности

p=i;

4. Запустить генетический алгоритм (ГА) для текущего ограничения по p;

5. Для каждой структуры, на каждой итерации ГА при вычислении «функции приспособленности» производим следующие действия:

1. Проверяем, была ли такая структура уже вычислена в store, если есть, извлекаем значения критерия и передаём в ГА;

2. Производим параметрическую идентификацию для текущей структуры-особи r, r. на основе критерия (2) и ограничения на параметры в виде алгебраического критерия устойчивости Рауса-Гурвица;

3. Добавляем модель в store;

4. Проверяем полученные параметры модели на устойчивость с помощью характеристического уравнения, и в случае устойчивости модели добавляем её в stab_store;

6. Сортируем stab_store по значению критерия от меньшего к большему и добавляем модель с меньшим критерием в выходной тезаурус store_p.

7. Увеличиваем границу поиска для ГА p=p+1;

8. Если при расширении пространства поиска p изменение критерия для оптимальной модели текущей p произойдет меньше чем на А по сравнению с предыдущим оптимальным значением т. е. p-1, то в этом случае останавливаем алгоритм;

9. Если пространство поиска p достигло максимально допустимой сложности, останавливаем алгоритм;

10. Переходим к п.4.

После остановки алгоритма выводим график зависимости критерия от p и получившийся тезаурус store_p, в котором модели уже отсортированы по сложности (размерности p) и указаны значения критерия (ошибки) для каждой модели.

На основе данного метода был разработан программный продукт «Структурно-параметрическая идентификация линейных динамических систем с помехами на входе и выходе», который позволяет фактически в автоматическом режиме производить структурно-параметрическую идентификацию ЛДС. В качестве среды для реализации данного продукта был выбран Matlab компании The MathWorks, одним из решающих преимуществ которого является богатое «математическое окружение» готовых функций и компонентов систем, которые необходимы для реализации данного метода.

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

Для решения задачи структурной идентификации в качестве основы был взят Genetic Algorithm

Рисунок 1. Упрощенная структура программы

Solver который реализовывает метод ГА и является частью пакета Matlab. В предоставленном виде данный пакет не адаптирован под решение задач целочисленного программирования (Integer programming). Для адаптации данного инструмента используется данный метод [13]:

1. Задается верхняя и нижняя целочисленная граница области, в которой будет происходить поиск.

2. Создается производная функция генерации начальной популяции, в которой гарантированно получаем целочисленное значение для каждой особи. Пример функции:

function Population = int_pop (GenomeLength, FitnessFcn, options) totalpopulation = sum (options. PopulationSize); range = options. PopInitRange; lower = range(1,:); span = range(2,:) - lower; Population = repmat (lower,totalpopulation,1) + round (repmat (span, totalpopulation, 1). *rand (totalpopulation, GenomeLength));

3. Создаём производную функцию мутации, в которое при мутировании получается строго целочисленное значение особи. Пример функции:

function mutationChildren = int_mutation (parents,options, GenomeLength, FitnessFcn, state,

thisScore, thisPopulation) shrink =.01; scale = 1; scale= scale - shrink * scale * state.Generation/ options.Generations; range = options. PopInitRange; lower=range(1,:); upper=range(2,:); scale = scale * (upper - lower); mutationPop = length(parents); mutationChildren = repmat (lower,mutationPop,1) + round(repmat(scale,mutationPop, 1).* rand (mutationPop, GenomeLength));

При данном подходе ГА реализованный в Matlab, будет производить пространственный поиск строго в рамках заданных границ (lower, upper), поиск будет проходить строго в целочисленной области.

Тестирование алгоритма.

В качестве системы для теста была выбранная следующая модель:

Число входных переменных x= 4, где запаздывание для каждого x будет: r1=3, r2=1, r3=1, r4=2, запаздывание по выходу r=2. Для нашего метода данная модель будет выглядеть в виде вектора r = [3,1,1,2,2], где каждое значение r соответствует порядку r. для каждого x, кроме последнего значения оно обозначает порядок г по выходу. Коэффициенты по выходу b = [0.8,-0.5,0.2] перенесены в правую часть модели (1) это необходимо для удобства реализации вычислений. Коэффициенты по входу a = [0.4,-0.5,-1,0.3,-0.2,0.6,0.4,-0.5,-1,0.3,-0.2] соответствующие сдвигу r. для каждого x. Сложность данной системы p=9, и для всех тестов общий объем выборки N = 10000.

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

• Метод наименьших квадратов;

• Рекуррентный метод инструментальных переменных [14], со следующим выбором вектора инструментальных переменных:

¥ i = yi-r-1 ••• У i-2r i w(—^ ...w(—2q

¡ w(d,í. ...w i i—rd

(d) i-2rd

• Разработанный нами критерий. Результат работы алгоритмов сравнивались по относительной погрешности моделирования:

&N =

z -100%,

где z = - вектор выходной ненаблюдае-

мой переменной,

т

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

2

z — z

Рисунок 2. График сходимости каждого критерия: а - разработанный критерий; б - инструментальные переменные; в - метод наименьших квадратов

Таблица 1. Сравнение критериев при р=9 (заданная сложность модели)

Так как результатом работы нашего метода является тезаурус моделей, возникает проблема сравнения моделей различных структур, поэтому для результирующей оценки будут браться модели одинаковой сложности р в каждом тезаурусе и модели при которых сработал критерий останова т.е Д<=1 [15, с. 356]. В качестве подтверждения данного высказывании провёлся тест на системе со следующими параметрами: о2 / о х = 0.2000; о, / о, = 0.2000;

ох 2 = [0.0396, 0.0403, 0.0415, 0.0396, 0.0397]

оу2 = 0.3497;

у = [0.1132, 0.1152, 0.1186, 0.1134, 0.1136].

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

Заключение

Таким образом, разработанный подход позволяет осуществлять структурно-параметрическую идентификацию ЛДС произвольной

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

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

1. Кацюба, О.А. Особенности применения МНК для оценивания линейных разностных операторов в задачах идентификации объектов управления / О.А. Кацюба, А.И. Жданов // Автоматика и телемеханика. - 1979. - №8. - С. 86-96.

°2/ О о,/О Погрешность 59,%

МНК Метод инструментальных переменных Разработанный критерий

0,2 0,2 8,61 8,39 7,32

0,2 0,5 12,8 13,44 7,37

0,2 0,7 15,3 16,1 7,56

0,5 0,2 18,7 18,69 10,72

0,7 0,2 28,87 29,19 7,63

0,5 0,5 23,13 24,09 10,21

0,5 0,7 25,11 22,89 10,35

0,7 0,5 33,15 33,26 12,42

0,7 0,7 31,31 31,65 10,76

Таблица 2. Сравнение критериев при Д < 1

°2/О о1 / О Погрешность 59 ,%

МНК Метод инструментальных переменных Разработанный критерий

0,2 0,2 3,83 3,98 1,17

0,2 0,5 4,40 4,62 2,73

0,2 0,7 5,06 7,23 3,84

0,5 0,2 17,5 17,45 2,82

0,7 0,2 28,24 28,07 3,67

0,5 0,5 18,51 18,31 3,13

0,5 0,7 18,97 18,37 7,51

0,7 0,5 29,58 29,88 4,38

0,7 0,7 31,57 32,53 5,93

2. Иванов, Д.В. Рекуррентная параметрическая идентификация многомерных линейных динамических систем с автокоррелированными помехами во входных и выходных сигналах / Д.В. Иванов, О.А. Кацюба // Вестник Самарского государственного технического университета, серия «Физико-математические науки науки». - 2011. - №4(25). - С. 102-109.

3. Иванов, Д.В. Рекуррентное оценивание билинейных динамических систем с помехами во входном и выходном сигналах / Д.В. Иванов, О.А. Усков // Известия Южного федерального университета, серия «Технические науки». - 2012. - №6 (131). - С. 187-192

4. Толчеек, В.О. Методы идентификации линейных одномерных динамических систем / В.О. Толчеек, Т.В. Яготкина. -М.: Московский энергетический институт, 1997. - 108 с.

5. Фатуев, В.А. Пакет программ структурной и параметрической идентификации линейных одномерных динамических систем / В.А. Фатуев, А.В. Каргин, В.М. Понятский // Труды второй всероссийской научной конференции «Проектирование инженерных научных приложений в среде Matlab». - М.: ИПУ РАН, 2004. - C. 715-762.

6. Фатуев, В.А. Структурно-параметрическая идентификация динамических систем : учебное пособие / В.А. Фатуев, А.В. Каргин, В.М. Поняский. - г.Тула: Издательство ТулГУ, 2003. - 156 с.

7. Структурно-параметрическая идентификация многомерных нестационарных динамических систем / В.А. Фатуев [и др.] // Труды III международной конференции «Идентификация систем и задач управления» (SICPR0'2004). - М.: ИПУ РАН, 2004. - С. 159-186.

8. Ljung, L. System Identification - Theory for the User / L. Ljung. - Prentice Hall, Upper Saddle River N.J., 2nd edition, 1999.

9. Анисимов, А.С. Идентификация порядка линейного разностного уравнения / А.С. Анисимов, В.Т. Кононов // Труды международной конференции «Идентификация систем и задачи управления» SICPRO «2000. Москва, 26-28 сентября 2000 г. Институт проблем управления им. В.А. Трапезникова РАН. - М: Институт проблем управления им. В.А. Трапезникова РАН, 2000. - 2534 с.

10. Кацюба, О.А. Теория идентификации стохастических динамических систем в условиях неопределённости: монография / О.А. Кацюба. - Самара: СамГУПС, 2008. - 119 с.

11. Емельянов, В. В. Теория и практика эволюционного моделирования / В.В. Емельянов, В.В. Курейчик, В.М. Курейчик. -М: Физматлит, 2003. - 432 с. - ISBN 5-9221-0337-7.

12. Matousek Jiri, Sharir Micha, Welzl Emo. A subexponential bound for linear programming, 1996, Algorithmica 16: 498-516, doi:10.1007/BF01940877.

13. A real coded genetic algorithm for solving integer and mixed integer optimization problems / Deep Kusum [etc.] // Applied Mathematics and Computation. - 2009. - 212(2). - P. 505-518.

14. Thil S. Contributions a l'identification de modeles avec des erreurs en les variables. Thesis of Doctorat, University of Henri Poincare, Nancy 1, 2007.

15. Льюнг, Л. Идентификация систем. Теория для пользователя / Л. Льюнг. - М.: Наука, 1991. - 432 с.

Сведение об авторе:

Энгельгардт Владислав Викторович, преподаватель кафедры мехатроника в автоматизированных производствах Самарского государственного университета путей сообщения (СамГУПС) 443066, г. Самара, 1-ый Безымянный переулок 18, e-mail: hexware@gmail.com