ON SOLUTION OF VARIATIONAL PROBLEM WITH GENETIC PROGRAMMING TECHNIQUES

详细

The authors offer to use genetic programming techniques to solve variational problem. Two solving variants of finding solution are described. Features of the problem solution are considered. Results of numerical experiments are given.

全文:

Вариационная задача состоит в нахождении функ- ции, обеспечивающей наибольшее (наименьшее) зна- чение заданного функционала [1]. При решении ва- риационной задачи обычно пользуются необходимым условием экстремума: с помощью уравнения Эйлера численно решают обыкновенное дифференциальное уравнение (ОДУ), в результате чего обычно получают экстремум функционала (если этот экстремум суще- ствует). Для решения ОДУ ставится краевая задача. В общем случае краевая задача не гарантирует сущест- вования и единственности решения, но даже если ре- шение существует и единственно, то далеко не каж- дое решение ОДУ может быть представлено в квадра- турах. При численном же решении краевой задачи Для решения вариационной задачи можно вос- пользоваться эволюционной теорией, а именно алго- ритмами генетического программирования (ГП) [2] как инструментом работы с бинарными деревьями. Данную задачу можно решить двумя методами: первый основан на использовании уравнения Эйлера, т. е. необходимого условия существования экстрему- ма, когда решение ОДУ ищется аналогично [3] (назо- вем его методом Эйлера), второй ищет функцию, дос- тавляющую наименьшее значение функционала, на- прямую (назовем его прямым методом). Решение вариационной задачи методом Эйлера. Этот метод сводит решение вариационной задачи к решению краевой задачи. При этом предполагается, возникают другие сложности: корректность постав- что функция F ( x, y( x), y′( x)) из (1) дифференцируема ленной задачи, определение устойчивости и сходимо- сти метода и т. д. В силу того что традиционные методы решения вариационной задачи имеют слабые стороны, необхо- димо разработать метод решения вариационной зада- чи, который мог бы объединить преимущества чис- ленного и аналитического решения, т. е. позволил бы численно получать функцию в символьном виде, ко- торая являлась бы экстремумом для заданного функ- трижды, а искомая функция дифференцируема дваж- ды [1]. Пусть из поставленной задачи получено обыкно- венное дифференциальное уравнение (уравнение Эй- лера) в виде F − d F = 0, y dx y ′ или ционала. Постановка задачи. Пусть задан функционал в Fy − Fxy′ − Fyy′ y′ − Fy′y′ y ′ = 0, и граничные условия (2) следующем виде: x1 y(x0) = Y0, y(xN) = YN, (3) v[ y( x)] = ∫ F ( x, y( x), y′( x))dx. (1) где y(x) – искомое решение; y′, y′′ – производные x0 Требуется определить функцию y(x), доставляющую минимум (максимум) этому функционалу. Предпола- гается, что на заданном интервале [x0; x1] функционал ограничен, а функция y(x) дважды непрерывно диф- ференцируема. Примечание 1. Если функционал задан в виде, отличном от выражения (1), то возможны два варианта решения зада- чи: преобразовать функционал, т. е. представить его в виде (1), или оставить его в первоначальном виде и изменить метод решения. Для простоты будем предполагать, что функционал может быть приведен к виду (1), хотя нет огра- ничений и для второго варианта решения задачи. функции y(x); Y0, YN – заданные вещественные кон- станты; N – объем выборки, т. е. количество точек разбиения отрезка. Требуется найти функцию y(x), удовлетворяющую уравнению (2) и граничным усло- виям (3). В теории обыкновенных дифференциальных урав- нений существование и единственность решения в общем случае не доказаны [4]. Но в поставленной задаче будем считать, что решение существует (для алгоритмов генетического программирования осталь- ные требования: единственности и непрерывности по входным данным – не являются необходимыми для решения задачи). *Работа выполнена при поддержке ФЦП «Исследования и разработки по приоритетным направлениям развития научно- технологического комплекса России на 2007–2013 годы» (НИР 2011-1.9-519.-005-042) и ФЦП «Научные и научно- педагогические кадры России» (НИР 2011-1.2.1-113-025). 19 Математика, механика, информатика Алгоритм генетического программирования дол- жен оперировать бинарными деревьями (в простей- шем варианте), представляющими собой потенциаль- ные решения вариационной задачи. Каждое бинарное дерево состоит из элементов функционального мно- жества (+, –, *, /, sin, cos, exp, log) и элементов терми- нального множества (x или компонентов вектора X, а также вещественных коэффициентов). В случае ре- шения специфических задач эти множества могут быть дополнены. В общем виде схему работы алгоритма можно представить следующим образом. Шаг 1. Инициализация начальной популяции (случайным образом). Шаг 2. Оценка пригодности каждого индивида. Шаг 3. Адаптация индивидов. Шаг 4. Применение генетических операторов к имеющимся индивидам для получения новой популя- ции. Шаг 5. Если критерий остановки не выполнен, то переход к шагу 2. На шаге 1 начальная популяция инициализируется случайным образом: до заданной начальной глубины функции пригодности от разницы между функцией и точками на границе. Примечание 2. В случае когда решение не единственное и существует конечное число решений, можно последова- тельно находить решения поставленной задачи, попутно штрафуя (например, добавив при вычислении пригодности дополнительное слагаемое) те решения, которые сходятся к уже найденным или известным решениям. Наличие тождества в уравнении (2) определяет функцию, являющуюся возможным решением задачи (1). Если функционал на этой функции достигает ми- нимума, то уравнение Эйлера превращается в тожде- ство. После проверки всех решений уравнения (2) (если решение не единственное) ищется локальный минимум (максимум) функционала (глобального экс- тремума функционал может и не достигать). Нужно учитывать и сложность полученной функ- ции, которая должна быть наименьшей: при идентич- ности графиков кривых предпочтение отдается функ- ции, которая имеет более короткую, лаконичную за- пись. Поэтому функцию пригодности можно предста- вить следующим образом: случайно выбираются элементы функционального множества и записываются в бинарное дерево, в ко- тором выбираются унарные операторы, и дальнейшее строительство дерева продолжается по одной ветви; Fit( Pk ) = 1 , 1 + E( Pk ) + K1 ⋅ D( Pk ) N ∑ Pk ( xi ) (4) на заданной глубине в дерево записываются случайно выбранные элементы терминального множества. Ко- личество индивидов в популяции также является па- E(Pk ) = i =1 + N (5) раметром работы алгоритма и должно в некоторой степени зависеть от поставленной задачи, так как слишком большая популяция требует значительных ресурсов для своей обработки, а слишком маленькая не дает точного и надежного результата. На шаге 2 происходит оценка каждого индивида согласно условиям поставленной задачи, при этом учитываются следующие факторы: – соответствие функции граничным условиям; – нарушение тождества при подстановке прове- ряемого решения задачи в уравнение (2); + K2 ⋅ ⎡⎣ y( x0 ) − Y0 + y( xN ) − YN ⎤⎦ , где Fit(Pk) – значение функции пригодности k-го ин- дивида; E(Pk) – ошибка аппроксимации, вычисляемая по всем точкам выборки, здесь N – объем выборки; K1 – коэффициент штрафа за сложность дерева; D(Pk) – число вершин дерева Pk; K2 – коэффициент штрафа за граничные условия; |Pk( xi )| – отклонение от нуля (например, среднеквадратическое) функции F из уравнения (2), получаемое при подстановке ре- шения Pk в ОДУ в точках (xi), выбранных из интер- вала равномерно, случайно или каким-либо специ- – сложность полученной функции. альным образом; Y ,Y – заданные граничные усло- Соответствие функции граничным условиям опре- 0 N деляет значительно более узкий круг функций, под- вия, y( x j ) – значение функции в точке xj. ходящих для решения вариационной задачи. При этом возможны три варианта: – через граничные точки не проходит ни одно ре- шение; – проходит единственное решение; – проходит несколько решений. В первом варианте задача не имеет решения при поставленных краевых условиях. Второй вариант яв- ляется наиболее желанным, когда решение существу- ет и единственно. В третьем варианте нужно найти несколько решений (этот случай будет рассмотрен ниже). Самым простым способом учета граничных условий является вычисление отклонения от них с некоторым коэффициентом. В более сложных слу- чаях может быть введена нелинейная зависимость При определении ошибки соответствия решения ОДУ требуется вычисление производных в точках выборки, которое можно выполнить двумя способа- ми: численным (с определенным порядком малости) [5] и аналитическим – путем построения дерева про- изводной функции и нахождения значения этой функции в точках выборки. Так как дифференцирова- ние (в отличие от интегрирования) всегда проводится по строгим правилам, то аналитически вычислить производную не представляет особого труда. При этом нетрудно заметить, что получение численной оценки производной требует меньше машинной памя- ти и времени, однако аналитически вычисленная про- изводная может быть точнее (при не слишком боль- шом количестве операторов в функции), поэтому ме- 20 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева тод вычисления производных должен выбираться с учетом специфики решаемой задачи. Примечание 3. В случае когда через граничные точки проходит не единственное решение, алгоритм может сой- тись к какому-то одному решению. Чтобы найти остальные решения, в функцию пригодности, т. е. в знаменатель вы- ражения (4), может быть добавлен дополнительный компо- нент, например такой: больше шансов стать родителями нового решения или быть отобранными в новую популяцию. Существуют разные виды селекции: пропорциональная, турнирная, ранговая и др. Скрещивание или клонирование (веро- ятности скрещивания и клонирования задаются) сле- дуют за селекцией: два отобранных индивида либо скрещиваются, либо один из них клонируется в новую популяцию. При скрещивании два дерева обменива- ются поддеревьями, отсеченными в случайно выбран- Pen( Pk ) = K3 ⋅ M 1 , N ∑(Sk ( xi ) − Fj ( xi )) ных точках. Это повышает вероятность того, что из деревьев с высокой пригодностью при скрещивании ∑ i =1 j =1 N где K3 – коэффициент штрафа за приближение к известному решению; М – количество частных решений (Fj), которые известны (или найдены на предыдущем этапе), и все они не являются искомым решением или являются не всеми иско- мыми; Sk – проверяемое решение. Примечание 4. В случае когда функционал представлен в виде, отличном от (1), и не сводится к нему из соображе- ний решения конкретной задачи, функция пригодности мо- жет быть представлена в виде, отличном от (4), (5), но обес- печивающем решение задачи. Этот случай требует отдель- ного исследования, поэтому в данной статье не рассматри- вается. На шаге 3 происходит адаптация индивидов, т. е. изменение частей бинарного дерева таким образом, чтобы значение функции пригодности было более высоким. Индивид, лучше адаптирующийся к услови- ям поставленной задачи, имеет лучшую пригодность, а следовательно, и более высокую вероятность быть отобранным для порождения потомков. С точки зре- ния теории ОДУ адаптация индивида означает изме- получатся деревья, превосходящие по пригодности родительские. Оператор мутации состоит в случайном изменении одного или нескольких узлов в дереве и рассматривается как средство восстановления генети- ческого разнообразия за счет случайного разброса решения по поисковому пространству. На шаге 5, если не выполнено условие остановки, следует переход на шаг 2 и повтор цикла. Условием остановки может быть как качественное ограничение – получение решения с заданной точностью, так и ко- личественное – ограниченность ресурсов (главным образом времени). Прямой метод решения вариационной задачи. При решении задачи (1) прямым методом схема алго- ритма генетического программирования остается прежней за исключением действий, связанных с вы- числением функции пригодности, поскольку пригод- ность оценивается по значению функционала. Если поставлена задача минимизации функциона- ла, то наилучшим решением является функция, кото- рая обеспечивает его минимум. Поэтому функцию пригодности можно представить следующим образом: нение его частей таким образом, чтобы решение име- ло меньшую погрешность, т. е. означает оптимизацию Fit( Pk ) = 1 , (6) 1 этого решения – изменение вещественных коэффици- ентов и/или функционального набора бинарного де- рева. Если зафиксировать все вершины дерева кроме какого-то одного узла, то его (дерево) можно рас- сматривать как функцию одной переменной. Теория оптимизации содержит много методов поиска экстре- мумов. Применительно к нашей задаче предлагается воспользоваться методом Хука–Дживса, отличаю- щимся простотой и эффективностью [6]. Бинарное дерево (решение ОДУ) в общем случае может содер- жать много вещественных коэффициентов и узлов- операций, и это делает целесообразным проведение лишь небольшого числа итераций. Для ускорения по- иска решения можно варьировать и функциональный набор, содержащийся в дереве, например при случай- ном изменении функционального терма, сохранить это изменение в случае улучшения или вернуться к старому терму в противном случае. На шаге 4 применяются генетические операторы селекции, рекомбинации (скрещивания), мутации. Оператор селекции – это оператор отбора индивидов с наибольшей пригодностью. Основная идея селекции состоит в том, что более пригодные индивиды имеют 1 + E( Pk ) + K1 ⋅ D( Pk ) + K2 ⋅ ∑ y( x j ) − Y j j =0 где E(Pk) – значение функционала на заданном интер- вале; остальные обозначения соответствуют форму- лам (4), (5). В случае если функционал представлен в виде (1), значение E(Pk) может быть определено численно, на- пример по формуле Симпсона на заданном множест- ве, когда последовательно находятся значения функ- ции и производных, а затем численное значение инте- грала. Полученное таким образом решение доставляет глобальный экстремум функционалу. Однако если экстремум не достижим, т. е. функ- ционал не ограничен снизу, на заданном множестве, то в этом случае можно говорить лишь о локальном экстремуме, найденном методом Эйлера. При этом апостериорное значение функционала ограничивается только машинной бесконечностью, что и определяет расходимость интеграла и, следовательно, неограни- ченность функционала. При прямом решении вариационной задачи в от- личие от решения с помощью уравнения Эйлера тре- буется больше машинных операций в силу того, что 21 Математика, механика, информатика сначала вычисляются значения функции и производ- ной (как и при решении уравнения Эйлера), а потом численно находится значение функционала, что зна- чительно увеличивает количество операций, необхо- димых для решения задачи, и время работы програм- мы. Однако в отличие от метода Эйлера прямой метод позволяет находить глобальный экстремум, поэтому далее мы будем рассматривать оба метода. Результаты численных экспериментов. Для ре- шения поставленной вариационной задачи была раз- работана программная система, обеспечивающая воз- можность решения как прямым методом, т. е. с по- мощью оценки значения функционала, так и с помо- щью оценки решения по уравнению Эйлера. При этом общая схема алгоритма отличается только вычисле- нием функции пригодности. В данной программной системе существует две группы настроек. Первая группа отвечает за текущую поставленную задачу и включает в себя собственно функционал (или ОДУ, представляющее уравнение Эйлера), граничные точки (интервал и граничные ус- ловия), функциональный набор, обеспечивающий представление решения. Вторая группа содержит на- стройки алгоритма генетического программирования: размер популяции и число поколений, ограничиваю- щие максимальное время работы алгоритма, вероят- ностные коэффициенты скрещивания, клонирования, мутации, вариант отбора (вид селекции). Оба набора настроек влияют на точность и время работы про- граммы, однако на надежность работы алгоритма наибольшее влияние оказывают настройки алгоритма ГП. Настройки задавались таким образом, чтобы со- блюдался некий компромисс между точностью, на- дежностью и временем работы. Для тестирования алгоритма были решены задачи, взятые из [1; 7] (табл. 1). Тестовые задачи были решены прямым методом и методом Эйлера. При решении каждым из методов для объективности оценки большинство настроек ал- горитма были неизменными, главным образом отли- чались только вычисленные значения функции при- годности. При этом для каждой задачи осуществлялся 20-кратный прогон алгоритма в силу его стохастично- сти, а результаты были усреднены. В случае когда функционал оказывался неограниченным (а это 5 из 10 задач), задача решалась только методом Эйлера. Пример полученного решения приведен в табл. 2. Там же отмечены номер поколения ГП, на котором было впервые найдено это решение, погрешность по- лучаемых решений и время работы алгоритма. Полученные решения можно разделить на три группы: – точные решения или решения, приводимые к точным элементарными преобразованиями без ок- ругления, например (((((x) · (x)) – (2,000)) · (x)) – (–1,000)); – условно точные решения, приводимые к точным элементарными преобразованиями с использованием округления, например (((sin((–4,000) · (x))) + ((exp(–34,800))/(19,100))) × × (–1,500)); – приближенные решения, т. е. решения, требую- щие более сложных преобразований и/или имеющие сложную громоздкую, не приводимую к точному ре- шению структуру дерева, например 61,9/(–39,1) · sin(log(31,900) + x). По результатам тестирования были вычислены средние показатели по всем задачам для каждой из групп при решении прямым методом и методом Эй- лера (табл. 3). Примеры решенных задач Таблица 1 № п/пФункция F ( x, y( x ), y′( x ))Граничные условияТочное решениеРешения, полученные алгоритмом 1(y′)2 + 12xy–3,0; 5,0 [–2; 2]x3–2x + 1((((x · x) – 2,0) · x) – (–1,0)) 2(y′)2 – y2–0,7; –0,28 [–3; 3]1,5sin(x) + 0,5cos(x)61,9/(–39,1) · sin(log(31,900) + x) 3(y′)2 + 2yy′ – 16y2–1,076; 1,076 [–1,7708 0,585398]1,5sin(4x))((sin((4,000) · (x))) · (1,500)) 4xy′ + (y′)2–4,0; –1,0 [–2; 4]–x2/4 + x – 1(((x · ((–4) + x)) + 4)/(–4,0)) 5(y′)2x2 + y′4,0; –0,667 [0, 4; 6]2/x – 1(((2,000) – (x))/(x)) 6(y′)2x + yy′–2,3; 1,79 [0, 1; 6]ln(x)(log(x)) 7(y′)2 – 2xy1,0; –1,0 [–3; 3](7x – x3)/6(((x · 45,5) + ((x · x) · (x · (–6,500))))/39,000) 8(y′)2 – y2 + 4y cos(x)0,141; 0,7 [–3; 3](x + 2)sin(x)log(exp(sin(x) · (2,0 + x))) 9(y′)2 – 2xy–0,81 2,43 [–2,5; 3,5]13x – x3/6 + 2((((x/(–6,0)) · x · x) + ((x + (2,0)) – (x/(–6,0)))) + x) 10(y′)2 + y2 + 2yex–4,19; –0,43[–1; 5]1,85914xex –0,5 + 0,5e–x((–4,1) – ((((x/((exp(–17,9)) · (–24,4) + cos(3,9))) – – exp(sin(17,4))) · exp(0,3))/exp(x))) + log(36,6) 22 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева Результат решения задачи с 20-кратным прогоном Таблица 2 Полученное решениеНомер поколения ГППогрешностьВремя, мин (((((x) · (x)) – (2,000)) · (x)) – (–1,000))1070,0952 ((1,000) + (((x) · ((x) · (x))) – ((x) + (x))))590,1143 ((1,000) – (((x) + (x)) – ((x) · ((x) · (x)))))480,1121 ((((1,000) + (((x) · (x)) · (x))) – (x)) – (x))1030,1153 ((((x) · (x)) · (x)) + ((1,000) – ((x) + (x))))780,1144 (((1,000) – (x)) + (((x) · ((x) · (x))) – (x)))370,1117 (((((x) · (x)) – (2,000)) · (x)) + (1,000))940,0980 (((x) · (–2,0)) + ((((x) · (x)) · (x)) + (cos((0,0) · (x)))))270,1420 (((1,000) + (((x) · (x)) · (x))) – ((x) + (x)))240,1113 ((1,000) + ((((x) · (x)) · (x)) – ((2,000) · (x))))320,1114 ((1,000) + ((x) · (((x) · (x)) – (2,000))))690,0929 (((((x) · (x)) · (x)) – ((x) + (–1,000))) – (x))370,1118 ((1,000) + ((((x) · (x)) – (2,000)) · (x)))290,0916 (((cos(–0,000)) – (x)) – ((x) – (((x) · (x)) · (x))))630,1250 ((1,000) – (((x) + (x)) – (((x) · (x)) · (x))))730,1148 ((((x) – ((3,000) · (x))) + (((x) · (x)) · (x))) – (–1,000))580,1332 ((((x) · (x)) · (x)) – (((x) + (x)) – (cos((x) – (x)))))270,1414 (((x) · ((x) · (x))) – ((–1,000) + ((x) + (x))))310,1117 (((1,000) + (((x) · (x)) · (x))) – ((x) + (x)))350,1123 ((((cos(0,000)) – (x)) + ((x) · ((x) · (x)))) – (x))440,1222 Среднее53,750,11131,3 Результаты анализа полученных решений Таблица 3 ПоказателиПрямой методМетод Эйлера Средний процент точных решений3562 Средний процент условно точных решений2525 Средний процент приближенных решений4013 Среднее количество поколений ГП12778 Среднее время работы, мин475321 Количество точных решений для наиболее сложной задачи010 Количество приближенных решений для наиболее сложной задачи10030 Анализируя полученные данные, можно сделать вывод о том, что и прямой метод, и метод Эйлера имеют право на существование. Прямой метод используется в случае поиска гло- бального экстремума, хотя в силу того что этот метод требует больше ресурсов и не отличается высокой надежностью (что можно объяснить накапливаемой вычислительной погрешностью при нахождении про- изводных, а потом при численном определении зна- чений функционала), его применение может быть обосновано только тем, что он не предполагает поис- ка производной более высокого порядка, как в методе Эйлера. Однако не каждую задачу можно решить та- ким методом, например в случае когда функционал неограничен сверху при поиске максимума (и снизу при поиске минимума), решения задачи не существу- ет (возможно, что это свойство данного метода может быть полезно для проверки ограниченности функцио- нала). Метод решения вариационной задачи с помощью уравнения Эйлера, т. е. перехода от вариационной задачи к краевой, показывает лучшие результаты и по времени работы, и по надежности, поэтому он пред- ставляется более полезным. Однако нужно отметить, что поставленная краевая задача не всегда имеет ре- шение, а если и имеет, то это решение может быть и не единственным или вообще не являться решением вариационной задачи. Поэтому такие задачи и их ре- шения должны быть подвергнуты дополнительному анализу. 23 Математика, механика, информатика Таким образом, предложен подход, основанный на алгоритме генетического программирования, позво- ляющий численно решать вариационную задачу в символьном виде, и разработана программа, реали- зующая этот подход, которая задает функционал в символьном виде и граничные условия и решает задачу с получением точной формулы, если она суще- ствует, или приближенного символьного выражения в случае когда решение нельзя выразить элементар- ными функциями. Возможно точное вычисление про- изводных путем построения производной функции или приближенное – с помощью разностных схем. Процесс решения может быть реализован двумя ме- тодами. Решения, полученные прямым методом, яв- ляются глобальными экстремумами, но в случае, ко- гда функционал неограничен, глобального экстрему- ма не существует и, следовательно, задачу решить нельзя. При решении методом Эйлера время работы алгоритма значительно сокращается, однако получе- ние решения является лишь необходимым условием существования минимума или максимума (возможно, локального). Поэтому предпочтение тому или иному методу решения должно отдаваться в зависимости от условий поставленной задачи, но в целом можно отметить, что метод Эйлера может использоваться чаще. Областью применения данного алгоритма может выступать любая отрасль науки и техники, так или иначе связанная с решением вариационной задачи. Однако этот алгоритм нужно рассматривать как до- полнительный инструмент в наборе уже имеющихся традиционных методов, таких как методы Ритца, Кан- торовича, Галеркина и др. Нужно также отметить, что традиционные методы ориентируются на решение известных типов уравнений, что и является ограниче- нием для их применения, поэтому эти методы и алго- ритм генетического программирования нельзя срав- нивать по точности или времени работы: очевидно, что традиционные методы будут быстрее, если ре- шаемая задача удовлетворяет их теоретическим по- ложениям. Тем не менее предложенный в данной ста- тье алгоритм заслуживает внимания в тех случаях, когда стандартные методы не могут дать желаемого результата.