Курсовая работа: Сравнительный анализ методов оптимизации

Министерство образования и науки Республики Казахстан

Карагандинский Государственный Технический Университет

Кафедра

Пояснительная записка

к курсовой работе по дисциплине: «Теория принятия решений»

Тема: Сравнительный анализ методов оптимизации

Выполнила

Студентка группы ______

______________________

Руководитель

______________________

Караганда-2009

Содержание

Введение

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

1 Прямые методы одномерной оптимизации

1.1 Метод дихотомии

1.2 Метод золотого сечения

2 Прямые методы безусловной оптимизации многомерной функции

2.1 Метод покоординатного циклического спуска

2.2 Метод Хука - Дживса

2.3 Метод правильного симплекса

2.4 Метод деформированного симплекса

3. Условная оптимизация

3.1 Метод преобразования целевой функции

3.2 Метод штрафных функций

4. Симплекс таблицы

Заключение

Список используемой литературы

Приложение А Листинг программ: Метод дихотомии, Метод золотого сечения, Метод покоординатного циклического спуска, Метод Хука – Дживса, Метод правильного симплекса

Приложение Б Листинг программы: Метод деформированного симплекса

Приложение В Листинг программы: Метод правильного трехмерного симплекса (максимизация объема фигуры)

Введение

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

Формулировка математической задачи оптимизации.

В достаточно общем виде математическую задачу оптимизации можно сформулировать следующим образом:

Минимизировать (максимизировать) целевую функцию с учетом ограничений на управляемые переменные.

Под минимизацией (максимизацией) функции n переменных f(x)=f(x1, ... ,xn) на заданном множестве U n-мерного векторного пространства En понимается определение хотя бы одной из точек минимума (максимума) этой функции на множестве U, а также, если это необходимо, и минимального (максимального) на U значения f(x).

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

f(x) -> min (max),

x принадлежит U,

где f(x) - целевая функция, а U - допустимое множество, заданное ограничениями на управляемые переменные.

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

Целью данного курсового проекта является изучение методов оптимизации функции. Методов одномерной оптимизации: метод дихотомии, золотого сечения; многомерной безусловной оптимизации: покоординатный циклический спуск, метод Хука – Дживса, правильный симплекс, деформированный симплекс, а также методов условной оптимизации Метод преобразования целевой функции, метод штрафных функций, табличный симплекс – метод.

1. Прямые методы одномерной оптимизации

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

f(x) -> min ,

x принадлежит [a, b].

Максимизация целевой функции эквивалента минимизации ( f(x) -> max ) эквивалентна минимизации противоположной величины ( -f(x) -> min ), поэтому, не умаляя общности можно рассматривать только задачи минимизации.

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

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

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

Рассмотрим наиболее распространенные на практике прямые методы поиска точки минимума. Самым слабым требованием на функцию f(x), позволяющим использовать эти методы, является ее унимодальность. Поэтому далее будем считать функцию f(x) унимодальной на отрезке [a, b].

Функция f(x) называется унимодальной на отрезке [a,b], если она непрерывна на [a,b] и существуют числа , , такие, что:

если , то на отрезке [a; ] функция f(x) монотонно убывает;

если , то на отрезке [; b] функция f(x) монотонно возрастает;

при

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

F(x)=8*cos(5*x)+ → min

a=2.7 b=3.9

И выбрано приближение ε=0,01, =0,02

1.1 Метод дихотомии

В этом методе точки x1 и х2 располагаются близко к середине очередного отрезка [а; b], т.е:

, (2.11)

где d > 0 – малое число. При этом отношение длин нового и исходного отрезков

близко к 1/2, этим и объясняется название метода.

Отметим, что для любых точек x1 и х2 величина t > 1/2, поэтому указанный выбор пробных точек объясняется стремлением обеспечить максимально возможное относительное уменьшение отрезка на каждой итерации поиска х*.

В конце вычислений по методу дихотомии в качестве приближенного значения х* берут середину последнего из найденных отрезков [а; b], убедившись предварительно, что достигнуто неравенство

.

Опишем алгоритм метода деления отрезка пополам.

Шаг 1. Определить x1 и х2 по формулам (2.11). Вычислить f (x1) и f (x2).

Шаг 2. Сравнить f (x1) и f (x2). Если , то перейти к отрезку [а; x2], положив b = x2 , иначе – к отрезку [x1; b], положив а = x1 .

Шаг 3. Найти достигнутую точность

Если , то перейти к следующей итерации, вернувшись к шагу 1. Если , то завершить поиск х*, перейдя к шагу 4.

Шаг 4. Положить

.

Замечания:

1. Число d из (2.11) выбирают на интервале (0;2e) с учетом следующих соображений:

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

б) при чрезмерно малом d сравнение значений f (x) в точках x1 и х2 , отличающихся на величину d, становится затруднительным. Поэтому выбор d должен быть согласован с точностью определения f (x) и с количеством верных десятичных знаков при задании аргумента х.

В таблице 1 приведено решение задания по варианту.

Таблица 1 - Метод дихотомии

|a-b|=0.001<= ε, x*=(a+b)/2=3.139282, f(x*)=-5.8074527

Рисунок 1 – Результат выполнения программы (Метод дихотомии).

1.2 Метод золотого сечения

Метод золотого сечения. Рассмотрим такое симметричное расположение точек x1 и х2 на отрезке [а; b], при котором одна из них становится пробной точкой и на новом отрезке, полученном после исключения части исходного отрезка. Использование таких точек позволяет на каждой итерации метода исключения отрезков, кроме первой, ограничиться определением только одного значения f (x), так как другое значение уже найдено на одной из предыдущих итераций.

Найдем точки x1 и х2 , обладающие указанным свойством.

Рассмотрим сначала отрезок [0; 1] и для определенности предположим, что при его уменьшении исключается правая часть этого отрезка. Пусть х2 = t, тогда симметрично расположенная точка х1 = 1–t (рис. 2.7).

Рис. 2.7. Определение пробных точек в методе золотого сечения

Пробная точка х1 отрезка [0; 1] перейдет в пробную точку х2¢ = 1–t нового отрезка [0; т]. Чтобы точки х2 = t, и х2¢ = 1–t делили отрезки [0; 1] и [0; t] в одном и том же отношении, должно выполняться равенство

или

,

откуда находим положительное значение

…

Таким образом,

х1 = 1–t = , .

Для произвольного отрезка [а; b] выражения для пробных точек примут вид

;

В таблице 2 приведено решение задания по варианту.

Опишем алгоритм метода золотого сечения.

Шаг 1. Найти х1 и х2 по формулам (2.15). Вычислить f (x1) и f (x2). Положить ,

.

Шаг 2. Проверка на окончание поиска: если en > e, то перейти к шагу 3, иначе – к шагу 4.

Шаг 3. Переход к новому отрезку и новым пробным точкам. Если f (x1) £ f (x2) то положить b=x2 , x2=x1 , f (x2) £ f (x1), x1=b–t(b–a) и вычислить f (x1), иначе – положить a=x1, x1= x2 , f (x1) = f (x2), x2=b+t(b–a) и вычислить f (x2). Положить en = ten и перейти к шагу 2.

Шаг 4. Окончание поиска: положить

, .

Результаты вычислений на остальных итерациях представлены в таблице 2 .

Таблица 2 - Метод золотого сечения

|a-b|=7.893370498E-3< ε, x*=(a+b)/2=3.1407091

f(x*)=-5.807126299

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

2 Прямые методы безусловной оптимизации многомерной функции

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

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

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

F(x1,x2)=a*x*y+(b*y+c*x)/x*y → min

a=5 b=3.5 c=2.5

x1=

x2=

2.1 Метод покоординатного циклического спуска

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

базисная точка переносится в

,

базисная точка переносится в

Циклы повторяются до тех пор, пока в ε окрестности найденной базисной точки будет улучшение функции. Решением поставленной задачи является точка в ε окрестности которой функция не принимает значение, лучшие, чем в этой точке.

Для решения поставленной задачи выбрано приближение ε=0,01 α=0,15

Таблица 3 - Метод покоординатного циклического спуска

α/2< ε, x1=0,2295 x2=0,4145 Z(x1,x2)= -0,8479571

2.2 Метод Хука – Дживса

Метод Хука и Дживса осуществляет два типа поиска - исследующий поиск и поиск по образцу. Первые две итерации процедуры показаны на рисунке 4.

Рисунок 4 – 1-поиск по образцу; 2- исследующий поиск вдоль координатных осей.

При заданном начальном векторе x1 исследующий поиск по координатным направлениям приводит в точку x2 . Последующий поиск по образцу в направлении x1- x2 приводит в точку y. Затем исследующий поиск, начинающийся из точки y, дает точку x3. Следующий этап поиска по образцу вдоль направления x3- x2 дает y*. Затем процесс повторяется.

Рассмотрим вариант метода, использующий одномерную минимизацию вдоль координатных направлений d1,..., dn и направлений поиска по образцу.

Начальный этап. Выбрать число eps > 0 для остановки алгоритма. Выбрать начальную точку x1, положить y1= x1, k=j=1 и перейти к основному этапу.

Основной этап.

Шаг 1. Вычилить lymj - оптимальное решение задачи минимизации f(yj+lym * dj) при условии lym принадлежит E1. Положить y[j+1]= yj+lymj*dj. Если j < n, то заменить j на j+1 и вернуться к шагу 1. Если j=n, то положить x[k+1] = y[n+1]. Если ||x[k+1] - xk|| < eps , то остановиться; в противном случае перейти к шагу 2.

Шаг 2. Положить d = x[k+1] - xk и найти lym - оптимальное решение задачи минимизации f(x[k+1]+lym*d) при условии lym принадлежит E1. Положить y1= x[k+1]+lym*d, j=1, заменить k на k+1 и перейти к шагу 1. Для решения поставленной задачи выбрано приближение ε=0,02, α=0,15

Таблица 4 - Метод Хука-Дживса

Т.к в ε окрестности полученной на 38 шаге точке мы не получаем улучшения (уменьшения значения) функции, то примем x1=0,426999999999999 x2=0,496999999999999,

Z(x1,x2)= -0,699011600000001.

2.3 Метод правильного симплекса

Правильный симплекс в пространстве En называется множество из n+1 равноудаленных друг от друга точек (вершин симплекса). В пространстве Е2 правильным симплексом является совокупность вершин равностороннего треугольника, Е3 – правильного тетраэдра.

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

Начальный этап. Выбрать параметр точности eps, базовую точку x0, ребро a и построить начальный симплекс. Вычислить f(x0).

Основной этап.

Шаг 1. Вычислить значения f(x) в вершинах симплекса x1,..., xn.

Шаг 2. Упорядочить вершины симплекса x0,..., xn так, чтобы f(x0)<=f(x1)<=...<=f(x[n-1])<=f(xn).

Шаг 3. Проверим на окончание поиска

,

где

Это одно из возможных условий останова. Его выполнении соответствует либо малому ребру a симплекса, либо попаданию точки минимума x* внутрь симплекса, либо тому и другому одновременно.

Если это условие выполнено, то вычисления прекратить, полагая x*= x0. В противном случае перейти к шагу 4.

Шаг 4. Найти xс и выполнить отражение вершины xn : y=2*xс- xn. Если f(y)<f(xn), то положить xn=y и перейти к шагу 2. Иначе - перейти к шагу 5.

Шаг 5. Перейти к новому правильному симплексу с вдвое меньшим ребром, считая базовой вершиной x0. Остальные n вершин симплекса найти по формуле xi=( xi+ x0)/2, i=1,...,n. Перейти к шагу 1.

Для решения поставленной задачи выбрано приближение ε=0,01, α=0,3

Таблица 5 - Метод симплекса

Т.к дальнейшее уменьшение α невозможно(α/2< ε) и в ε окрестности полученной на 31 шаге точке мы не получаем улучшения (уменьшения значения) функции, то примем x=0,248249999999998 и y=0,408289729858682 Z(x,y)= -0,847811621576176.

2.4 Метод деформированного симплекса

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

z1= xc - a( xc - xn), 0<a<1;

z2 = xc + a( xc - xn), 0<a<1;

z3 = xc + b( xc - xn), b приближенно равно 1;

z4 = xc + g( xc - xn), g<1.

Геометрическая иллюстрация этих процедур для двумерного пространства приведена на рисунке 7.

Новые симплексы полученные в результате процедуры сжатия (a,b); отражения (c); растяжения (d)

Так как величина a принадлежит интервалу (0;1), то выбор точек z1 и z2 соответствует сжатию симплекса; b приближенно равно 1, поэтому выбор точки z3 соответствует отражению, а g>1 и выбор точки z4 приводит к растяжению симплекса.

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

Алгоритм метода поиска точки минимума функции по деформируемому симплексу

Начальный этап. Выбрать параметр точности eps, параметры a, b и g, базовую точку x0 , параметр a и построить начальный симплекс. Вычислить значение функции f(x0).

Основной этап. Шаг 1. Вычислить значения функции в вершинах симплекса x1,..., xn.

Шаг 2. Упорядочить вершины симплекса x0,..., xn так, чтобы f(x0)<=f(x1)<=...<=f(x[n-1])<=f(xn).

Шаг 3. Проверить условие (1/n)Sum[f(xi)-f(x0)]^2 < e^2, i=[1,n].

Это одно из возможных условий останова. При его выполнении "дисперсия" значений f(x) в вершинах симплекса становится меньше e2, что, как правило, соответствует либо малому ребру a симплекса, либо попаданию точки минимума x* внутрь симплекса, либо тому и другому одновременно.Если это условие выполнено, то вычисления прекратить, полагая x*= x0. В противном случае перейти к шагу 4.

Шаг 4. Найти xс и пробные точки zk , k=1,...,4 по формулам (2). Найти f(z*)=minf(zk). Если f(z*)<f(xn), то положить xn = z* и перейти к шагу 2. Иначе - перейти к шагу 5.

Шаг 5. Уменьшить симплекс, полагая xi=( xi+ x0)/2, i=1,...,n перейти к шагу 1.

На практике хорошо зарекомендовал себя следующий набор параметров a=1/2, b=1, g=2, поэтому он и был использован в программе.

Таблица 5 – Метод деформированного симплекса

T.к.

< ε,

то примем x=0,237037034153931 и y=0,407407409218273 Z(x,y)= -0,848148148148148

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

3. Условная оптимизация

Задача условной оптимизации в общем случае записывается в известном виде:

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

На рисунке 12 представлена фигура, объем, которой необходимо максимизировать при заданной площади поверхности

Рисунок 12 – Фигура для максимизации объема при заданной площади поверхности

Найдем полную площадь поверхности данной фигуры(без верхней поверхности):

,

найдем объем фигуры:

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

Эту задачу можно решить двумя методами:

Метод преобразования целевой функции,

метод штрафных функций.

3.1 Метод преобразования целевой функции

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

V = 4/3∙a2∙h2+7/3∙h1∙a2 → max (1)

S = 6∙a∙h1+4∙h2∙a (2)

Выразим a из (2) и подставим в (1), получим:

V = s2∙(4∙h2+7∙h1)/3∙(6∙h1+4∙h2)2

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

Возьмем, например, метод правильного симплекса, и зададим начальные условия: а=1м, h1=3м, h2=4м, s=34м. Для метода симплекса выберем точность ε=0,001.

Т.е максимальный объем V=12,7151461307724, при заданной площади получается при h1 = 2,946875, и h2 = 3,83229490168751

3.2 Метод штрафных функций

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

f(x) -> min;

gi(x) 0, i 1, ..., k;

hj(x) 0, j 1, ..., m;

a x b.

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

F(x,a) f(x) +rS(x)

Здесь f(x) - целевая функция задачи оптимизации; S(x) - специальным образом выбранная функция штрафа,где r— множитель, значения которого можно изменять в процессе оптимизации.. Точку безусловного минимума функции F(x, a) будем обозначать через х(а).

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

3.2 Методы штрафных функций

Эти методы применяются для решения задач нелинейного программирования с ограничениями-неравенствами.

В рассматриваемых методах функции Ф(x, а) подбирают такими, чтобы их значения неограниченно возрастали при приближении к границе допустимой области G (Рисунок 14). Иными словами, приближение к границе “штрафуется” резким увеличением значения функции F(x, а). На границе G построен “барьер”, препятствующий нарушению ограничении в процессе безусловной минимизации F(x, a). Поиск минимума вспомогательной функции F(x, а) необходимо начинать с внутренней точки области G .

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

Здесь dG -граница области G.

Рисунок 14 - Внутренняя штрафная функция

Методы внешних штрафных функций

Данные методы применяются для решения задачи оптимизации при наличии как ограничений-неравенств, так и ограничений-равенств.

В рассматриваемых методах функции Ф(х, а) выбирают такими, что их значения равны нулю внутри и на границе допустимой области G, а вне ее -положительны и возрастают тем больше, чем сильнее нарушаются ограничения (Рисунок 15). Таким образом, здесь “штрафуется” удаление от допустимой области G.

Рисунок – 15 Внешняя штрафная функция

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

Для данного курсового проекта штрафная функция для объема данной фигуры имеет вид:

,

где - параметр штрафа, С – полная площадь поверхности, заданная изначально, V(a,h1,h2) = 4/3∙a2∙h2+7/3∙h1∙a2, S(a,h1,h2) = 6∙a∙h1+4∙h2∙a.

Задача была решена методом правильного трехмерного симплекса.

Мы видим, что при увеличении значения параметра штрафа, значение функции уменьшается (ухудшается), а при уменьшении – увеличивается (улучшается).

4. Симплекс таблицы

Для его применения необходимо, чтобы знаки в ограничениях были вида «меньше либо равно», а компоненты вектора b - положительны.

Алгоритм решения сводится к следующему:

Приведение системы ограничений к каноническому виду путём введения дополнительных переменных для приведения неравенств к равенствам.

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

Формируется симплекс-таблица.

Рассчитываются симплекс-разности.

Принимается решение об окончании либо продолжении счёта.

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

7 На каждой итерации определяется вектор, вводимый в базис, и вектор, выводимый из базиса. Пересчитывается таблица.

Дана функция вида:

f(x)=4x1+2x2

Подберем k геометрическим способом решения так, чтобы область допустимых значений была пятиугольником. k=7

Рисунок – 16 Область допустимых значений

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

4у1+2у2+0у3 +0у4 +0у5

=(0;0;8;12;7) – начальные допустимые базисные решения

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

Вводим в базис у1 , а выводим из базиса у4.

Вводим в базис у2 , а выводим из базиса у5.

Т.к. f<0, то останавливаемся на второй итерации.

Исходя из графика ОДЗ, можно определить, что оптимальным решением является отрезок прямой , входящий в

ОДЗ, проверим: 2,5*2+2=7.

x1 = 2,5, x2 = 2 f(x)=14.

Заключение

Целью данного курсового проекта было изучение методов оптимизации функции. Методов одномерной оптимизации: метод дихотомии, золотого сечения; многомерной безусловной оптимизации: покоординатный циклический спуск, метод Хука – Дживса, правильный симплекс, деформированный симплекс, а также методов условной оптимизации Метод преобразования целевой функции, метод штрафных функций, табличный симплекс – метод.

Список используемой литературы

1. А.Г.Трифонов. Постановка задачи оптимизации и численные методы ее решения;

2. Б. Банди. Методы оптимизации. Вводный курс., 1988;

3. Мендикенов К.К. Лекции

Приложение А

using System;

using System.Collections.Generic;

using System.ComponentModel;

using System.Data;

using System.Drawing;

using System.Text;

using System.Windows.Forms;

namespace lab1

{

public partial class Form1 : Form

{

class global

{

private global() { }

public static double a = 0.64;

public static double b = 1.77;

public static double e = 0.0001;

public static double al = 0.00001;

public static double x = 0;

public static double y = 0;

public static int iter = 0;

}

public Form1()

{ InitializeComponent(); }

private void textBox1_TextChanged(object sender, EventArgs e)

{global.e = Convert.ToDouble(textBox1.Text); }

private void textBox2_TextChanged(object sender, EventArgs e)

{ global.al = Convert.ToDouble(textBox2.Text); }

public double F(double x)

{ return (Math.Pow((2.5 - x), 2) + 3.1 * x); }

public double Z(double x, double y)

{ return (2.5 * Math.Pow(x, 2) + 2 * x * y + 3.1 * Math.Pow(y, 2) -2* x-3*y); }

public double Dixotom()

{

global.iter = 1;

global.a = Convert.ToDouble(textBox4.Text);

global.b = Convert.ToDouble(textBox3.Text);

richTextBox1.Text = richTextBox1.Text+"a="+Convert.ToString(global.a)+"; b="+Convert.ToString(global.b)+(char)13;

while (true)

{

double x1 = (global.a+global.b)/2-global.al;

double x2 = (global.a + global.b) / 2 + global.al;

if (F(x1) < F(x2)) global.b = x2; else global.a = x1;

richTextBox1.Text = richTextBox1.Text + Convert.ToString(global.iter) + ") x1=" + Convert.ToString(x1) + "; x2=" + Convert.ToString(x2) + "; f(x1)=" + Convert.ToString(F(x1)) + "; f(x2)=" + Convert.ToString(F(x2)) + "; a=" + Convert.ToString(global.a) + "; b=" + Convert.ToString(global.b) + (char)13;

global.iter++;

if (Math.Abs(global.b - global.a) < global.e) break;

} return (global.a + global.b) / 2;

}

public double Zolot()

{ global.iter = 1;

global.a = Convert.ToDouble(textBox4.Text);

global.b = Convert.ToDouble(textBox3.Text);

richTextBox1.Text = richTextBox1.Text + "a=" + Convert.ToString(global.a) + "; b=" + Convert.ToString(global.b) + (char)13;

double x2 = global.a+0.618*(global.b - global.a) ;

double x1 = global.a + (1-0.618) * (global.b - global.a);

while (true)

{

if (Math.Abs(global.b - global.a) < global.e) break;

richTextBox1.Text = richTextBox1.Text + Convert.ToString(global.iter) + ") a=" + Convert.ToString(global.a) + "; b=" + Convert.ToString(global.b) + "; x1=" + Convert.ToString(x1) + "; x2=" + Convert.ToString(x2) + "; f(x1)=" + Convert.ToString(F(x1)) + "; f(x2)=" + Convert.ToString(F(x2)) + (char)13;

if (F(x2) > F(x1))

{ global.b = x2; x2 = x1; x1 = global.a + 0.372 * (global.b - global.a); }

else { global.a = x1; x1 = x2; x2 = global.a + 0.618 * (global.b - global.a); }

global.iter++;

}

return (global.a + global.b) / 2;

}

private void button1_Click(object sender, EventArgs e)

{ richTextBox1.Text = "";

global.al = Convert.ToDouble(textBox2.Text);

global.e = Convert.ToDouble(textBox1.Text);

if (radioButton1.Checked) global.x = Dixotom();

if (radioButton2.Checked) global.x = Zolot();

label2.Text = "Минимум: x*=" + Convert.ToString(global.x) + "; y(x*)=" + Convert.ToString(F(global.x)) + ", число итераций: "+Convert.ToString(global.iter-1);

}

public void Spusk(double x,double y)

{

while (true)//идем вправо

{ x = x + global.al; if (Z(x, y) > Z(x - global.al, y)) break; global.iter++;

richTextBox2.Text = richTextBox2.Text + Convert.ToString(global.iter)+ ") x=" + Convert.ToString(x) + "; y=" + Convert.ToString(y) + "; z(x,y)=" + Convert.ToString(Z(x, y)) + "; al=" + Convert.ToString(global.al) + (char)13;

x = x - global.al;//возвращаемся на неудачный шаг

while (true)//идем влево

{ x = x - global.al; if (Z(x, y) > Z(x + global.al, y)) break; global.iter++;

richTextBox2.Text = richTextBox2.Text + Convert.ToString(global.iter) + ") x=" + Convert.ToString(x) + "; y=" + Convert.ToString(y) + "; z(x,y)=" + Convert.ToString(Z(x, y)) + "; al=" + Convert.ToString(global.al) + (char)13; }

x = x + global.al;//возвращаемся на неудачный шаг

global.x=x; global.y=y;

SpuskV(x, y);

}

public void SpuskV(double x, double y)

{

while (true)//идем вверх

{ y = y + global.al; if (Z(x, y) > Z(x, y - global.al)) break; global.iter++;

richTextBox2.Text = richTextBox2.Text + Convert.ToString(global.iter) + ") x=" + Convert.ToString(x) + "; y=" + Convert.ToString(y) + "; z(x,y)=" + Convert.ToString(Z(x, y)) + "; al=" + Convert.ToString(global.al) + (char)13; }

y = y - global.al;//возвращаемся на неудачный шаг

while (true)//идем вниз

{ y = y - global.al; if (Z(x, y) > Z(x, y + global.al)) break; global.iter++;

richTextBox2.Text = richTextBox2.Text + Convert.ToString(global.iter) + ") x=" + Convert.ToString(x) + "; y=" + Convert.ToString(y) + "; z(x,y)=" + Convert.ToString(Z(x, y)) + "; al=" + Convert.ToString(global.al) + (char)13; }

y = y + global.al;//возвращаемся на неудачный шаг

global.x = x; global.y = y;

if (global.al/2 > global.e) { global.al = global.al / 2; Spusk(x, y); }

}

public void Hyg(double x, double y)

{ while (true)

{int min=Vibor(x, y);

if (min == 1) { x = x + 2 * global.e; y = y + 2 * global.e; if (Z(x - 2 * global.e, y - 2 * global.e) < Z(x, y)) break; }

if (min == 2) { x = x - 2 * global.e; y = y + 2 * global.e; if (Z(x + 2 * global.e, y - 2 * global.e) < Z(x, y)) break; }

if (min == 3) { x = x - 2 * global.e; y = y - 2 * global.e; if (Z(x + 2 * global.e, y + 2 * global.e) < Z(x, y)) break; }

if (min == 4) { x = x + 2 * global.e; y = y - 2 * global.e; if (Z(x - 2 * global.e, y + 2 * global.e) < Z(x, y)) break; }

global.iter++;

richTextBox2.Text = richTextBox2.Text + Convert.ToString(global.iter) + ") x=" + Convert.ToString(x) + "; y=" + Convert.ToString(y) + "; z(x,y)=" + Convert.ToString(Z(x, y)) +(char)13; }

global.x = x; global.y = y;

}

public int Vibor(double x, double y)

{ int min = 0;

if (Z(x + global.e, y + global.e) < Z(x, y)) min = 1;

if (Z(x + global.e, y - global.e) < Z(x, y)) min = 2;

if (Z(x - global.e, y - global.e) < Z(x, y)) min = 3;

if (Z(x - global.e, y + global.e) < Z(x, y)) min = 4;

return min; }

public void Sym(double x, double y)

{ double x0 = x; double y0 = y; double x1 = x0 + global.al; double y1 = y0;

double x2 = x0 + (global.al) / 2; double y2 = y0 + global.al * Math.Sin(60);

richTextBox2.Text = richTextBox2.Text + Convert.ToString(global.iter) + ") z(x0,y0)=" + Convert.ToString(Z(x0, y0)) + " z(x1,y1)=" + Convert.ToString(Z(x1, y1)) + " z(x2,y2)=" + Convert.ToString(Z(x2, y2)) + " al=" + Convert.ToString(global.al) + (char)13;

while (true)

{ //поиск наименьшего

double mx0 = x0; double my0 = y0; double mx1 = x1; double my1 = y1; double mx2 = x2; double my2 = y2;

double z1 = Z(mx0, my0); double z2 = Z(mx1, my1); double z3 = Z(mx2, my2);

if ((z1 < z2) && (z2 < z3) && (z3 > z1)) { x0 = mx0; x1 = mx1; x2 = mx2; y0 = my0; y1 = my1; y2 = my2; }

if ((z1 < z2) && (z2 > z3) && (z3 > z1)) { x0 = mx0; x1 = mx2; x2 = mx1; y0 = my0; y1 = my2; y2 = my1; }

if ((z1 > z2) && (z2 < z3) && (z3 < z1)) { x0 = mx1; x1 = mx2; x2 = mx0; y0 = my1; y1 = my2; y2 = my0; }

if ((z1 > z2) && (z2 < z3) && (z3 > z1)) { x0 = mx1; x1 = mx0; x2 = mx2; y0 = my1; y1 = my0; y2 = my2; }

if ((z1 < z2) && (z2 > z3) && (z3 < z1)) { x0 = mx2; x1 = mx0; x2 = mx1; y0 = my2; y1 = my0; y2 = my1; }

if ((z1 > z2) && (z2 > z3) && (z3 < z1)) { x0 = mx2; x1 = mx1; x2 = mx0; y0 = my2; y1 = my1; y2 = my0; }

//проверка на выход

if (global.al <= global.e) break;

while (true)

{ //отражение относительно 3

double kx= (x0+x1)-x2; double ky = (y0 + y1) - y2;

if (Z(x2, y2) > Z(kx, ky)) { x2 = kx; y2 = ky; global.iter++; break; }

//отражение относительно 2

kx = (x0 + x2) - x1; ky = (y0 + y2) - y1;

if (Z(x1, y1) > Z(kx, ky)) { x1 = kx; y1 = ky; global.iter++; break; }

//отражение относительно 1

kx = (x1 + x2) - x0; ky = (y1 + y2) - y0;

if (Z(x0, y0) > Z(kx, ky)) { x0 = kx; y0 = ky; global.iter++; break; }

//уменьшаем треугольник

global.al = global.al / 2;

x1 = (x0 + x1) / 2; y1 = (y0 + y1) / 2;

x2 = (x0 + x2) / 2; y2 = (y0 + y2) / 2;

}

richTextBox2.Text = richTextBox2.Text + Convert.ToString(global.iter) + ") x0=" + Convert.ToString(x0) + " x1=" + Convert.ToString(x1) + " x2=" + Convert.ToString(x2) + "; y0=" + Convert.ToString(y0) + " y1=" + Convert.ToString(y1) + " y2=" + Convert.ToString(y2) + " z(x0,y0)=" + Convert.ToString(Z(x0, y0)) + " z(x1,y1)=" + Convert.ToString(Z(x1, y1)) + " z(x2,y2)=" + Convert.ToString(Z(x2, y2)) + " al=" + Convert.ToString(global.al) + (char)13; }

global.x = x0; global.y = y0; }

private void button2_Click(object sender, EventArgs e)

{ global.iter = 0;

global.al = Convert.ToDouble(textBox7.Text);

global.e = Convert.ToDouble(textBox8.Text);

if (radioButton4.Checked) {Spusk(Convert.ToDouble(textBox6.Text),Convert.ToDouble(textBox5.Text)); }

if (radioButton3.Checked)

Hyg(Convert.ToDouble(textBox6.Text), Convert.ToDouble(textBox5.Text));

if (radioButton5.Checked) Sym(Convert.ToDouble(textBox6.Text), Convert.ToDouble(textBox5.Text));

label9.Text = "Минимум: (" + Convert.ToString(global.x) + ";" + Convert.ToString(global.y) + " f(x,y)=" + Convert.ToString(Z(global.x, global.y)) + "), число итераций: " + Convert.ToString(global.iter); } }}

Приложение Б

procedure TForm1.Button1Click(Sender: TObject);

begin

a:=false;

l:=0.05;

al:=1; e:=0.01; gm:=2; bt:=0.5;

x:=1.16166; y:= 1.15185;

iter:=0;

xl:= x; yl:= y; xg:= xl + l; yg:= yl;

xh:= xl + l / 2; yh:= yl + l * Sin(60);

Sym(x,y);

ZZ(x,y); //считаем значение функции в найденной точке

end;

procedure TForm1.Sym(x,y:real);

label

ext,B;

begin

//поиск наименьшего

B: mx0:= xl; my0:= yl; mx1:= xg; my1:= yg; mx2:= xh; my2:= yh;

ZZ(mx0, my0); z1:=z; ZZ(mx1, my1); z2:=z; ZZ(mx2, my2); z3:=z;

if ((z1 < z2) and (z2 < z3) and (z3 > z1)) then begin xl:= mx0; xg:= mx1; xh:= mx2; yl:= my0; yg:= my1; yh:= my2; end;

if ((z1 < z2) and (z2 > z3) and (z3 > z1)) then begin xl:= mx0; xg:= mx2; xh:= mx1; yl:= my0; yg:= my2; yh:= my1; end;

if ((z1 > z2) and (z2 < z3) and (z3 < z1)) then begin xl:= mx1; xg:= mx2; xh:= mx0; yl:= my1; yg:= my2; yh:= my0; end;

if ((z1 > z2) and (z2 < z3) and (z3 > z1)) then begin xl:= mx1; xg:= mx0; xh:= mx2; yl:= my1; yg:= my0; yh:= my2; end;

if ((z1 < z2) and (z2 > z3) and (z3 < z1)) then begin xl:= mx2; xg:= mx0; xh:= mx1; yl:= my2; yg:= my0; yh:= my1; end;

if ((z1 > z2) and (z2 > z3) and (z3 < z1)) then begin xl:= mx2; xg:= mx1; xh:= mx0; yl:= my2; yg:= my1; yh:= my0; end;

Richedit1.Lines.Add(IntToStr(iter));

Richedit1.Lines.Add(FloatToStr(xl)+' '+FloatToStr(yl)+' '+FloatToStr(zl));

Richedit1.Lines.Add(FloatToStr(xg)+' '+FloatToStr(yg)+' '+FloatToStr(zg));

Richedit1.Lines.Add(FloatToStr(xh)+' '+FloatToStr(yh)+' '+FloatToStr(zh));

Richedit1.Lines.Add('');

x0:=(xl+xg)/2; y0:=(yl+yg)/2;

xr:=(1+al)*x0-al*xh; yr:=(1+al)*y0-al*yh;

ZZ(xl,yl); zl:=z; ZZ(xg,yg); zg:=Z; ZZ(xh,yh); zh:=z; ZZ(xr,yr); zr:=z;

if zr<zl then

{1} begin

//растяжение

xe:=gm*xr+(1-gm)*x0; ye:=gm*yr+(1-gm)*y0; ZZ(xe,ye); ze:= z;

if ze<zl then

{2} begin

xh:=xe; yh:=ye; ZZ(xh,yh); zh:=z;

xl:=gm*xl+(1-gm)*x0; yl:=gm*yl+(1-gm)*y0;ZZ(xl,yl); zl:=z;

xg:=gm*xg+(1-gm)*x0; yg:=gm*yg+(1-gm)*y0; ZZ(xg,yg); zg:=z;

Shod(xl,xg,xh,yl,yg,yh);

if a=true then goto ext else begin inc(iter); goto B; end;

{2} end

else

{3} begin

xh:=xr; yh:=yr; ZZ(xh,yh); zh:=z;

Shod(xl,xg,xh,yl,yg,yh);

if a=true then goto ext else begin inc(iter); goto B; end;

{3} end;

{1} end;

if ((zr>zl)and(zr<=zg)) then

begin

xh:=xr; yh:=yr; ZZ(xh,yh); zh:=z;

Shod(xl,xg,xh,yl,yg,yh);

if a=true then goto ext else begin inc(iter); goto B; end;

end;

if ((zr>zl)and(zr>zg)) then

{4E} begin

if zr<zh then begin xh:=xr; yh:=yr; ZZ(xh,yh); zh:=z; end;

//сжатие

xc:=bt*xh+(1-bt)*x0; yc:=bt*yh+(1-bt)*y0; ZZ(xc,yc); zc:=z;

if zc<xh then

{5Ж} begin

xh:=xc; yh:=yc; ZZ(xh,yh); zh:=z;

xl:=bt*xl+(1-bt)*x0; yl:=bt*yl+(1-bt)*y0;ZZ(xl,yl); zl:=z;

xg:=bt*xg+(1-bt)*x0; yg:=bt*yg+(1-bt)*y0; ZZ(xg,yg); zg:=z;

Shod(xl,xg,xh,yl,yg,yh);

if a=true then goto ext else begin inc(iter); goto B; end;

{5} end

{6З} else begin

//уменьшение

l:=l/2;

xh:=(xl+xh)/2; xh:=(xl+xh)/2; ZZ(xh,yh); zh:=z;

xg:=(xl+xg)/2; yg:=(yl+yg)/2; ZZ(xg,yg); zg:=z;

{6} end;

Shod(xl,xg,xh,yl,yg,yh);

if a=true then goto ext else begin inc(iter); goto B; end;

{4} end;

ext: Richedit1.Lines.Add(FloatToStr(xl));

Richedit1.Lines.Add(FloatToStr(yl));

Richedit1.Lines.Add(FloatToStr(zl));

Richedit1.Lines.Add(IntToStr(iter));

end;

procedure TForm1.ZZ(x,y:real);

begin

z:= 2.5 * x*x + 2 * x * y + 3.1 * y*y -2* x-3*y;

end;

procedure TForm1.Shod(xl,xg,xh,yl,yg,yh:real);

var

sigma:real;

begin

ZZ(xl,yl); zl:=z; ZZ(xg,yg); zg:=Z; ZZ(xh,yh); zh:=z;

sigma:=sqrt((sqr(zl-((zl+zg+zh)/2+1))+sqr(zl-((zl+zg+zh)/2+1))+sqr(zg-((zl+zg+zh)/2+1))+sqr(zh-((zl+zg+zh)/2+1)))/3); if sigma<e then a:=true;end;

Приложение В

unit Unit1;

type

procedure Button1Click(Sender: TObject);

procedure ZZ(x,y,z:real);

procedure Sym(x,y,z:real);

end;

var

X, y,z:real;

z_:real; //значение функции

iter:integer; //число итераций

s,gm,x0,x1,x2,x3:real; //координаты точек

y0,y1,y2,y3:real;

z0,z1,z2,z3:real; // треугольника

al, e:real; // длина стороны симпликса(треугольника)

kx,ky,kz, zz1,zz2:real; //координаты точки k

a,b:boolean; //для цикла

implementation

procedure TForm1.Sym(x,y,z:real);

label

l;

var

a:array[1..4,1..4] of real;//z()xyz ; 1234

k,i:integer;

buf:real;

changed:boolean;

{1} begin

x0:= x; y0:= y;z0:=z;

x1:= x0 + al; y1:= y0;z1:=z;

x2:= x0 + al / 2; y2:= y0 + al * Sin(60);z2:=z;

x3:= x0 + al/ 2; y3:=y0 ;z3:=z0 + al * Sin(60) ;

while (true) do

//for i:=1 to 100 do

{2} begin

ZZ(x0, y0,z0); a[1,1]:=z_; a[2,1]:=x0; a[3,1]:=y0; a[4,1]:=z0;

ZZ(x1, y1,z1); a[1,2]:=z_; a[2,2]:=x1; a[3,2]:=y1; a[4,2]:=z1;

ZZ(x2, y2,z2); a[1,3]:=z_; a[2,3]:=x2; a[3,3]:=y2; a[4,3]:=z2;

ZZ(x3, y3,z3); a[1,4]:=z_; a[2,4]:=x3; a[3,4]:=y3; a[4,4]:=z3;

//сортировка

repeat

changed:=FALSE;

for k:=1 to 3 do

if a[1,k] > a[1,k+1] then

begin

buf := a[1,k]; a[1,k] := a[1,k+1]; a[1,k+1] := buf;

buf := a[2,k]; a[2,k] := a[2,k+1]; a[2,k+1] := buf;

buf := a[3,k]; a[3,k] := a[3,k+1]; a[3,k+1] := buf;

buf := a[4,k]; a[4,k] := a[4,k+1]; a[4,k+1] := buf;

changed := TRUE;

end;

until not changed;

x0:=a[2,1]; y0:=a[3,1]; z0:=a[4,1];

x1 :=a[2,2]; y1:=a[3,2]; z1:=a[4,2];

x2:=a[2,3]; y2:=a[3,3]; z2:=a[4,3];

x3 :=a[2,4]; y3:=a[3,4]; z3:=a[4,4];

ZZ(x3, y3,z3);

//проверка на выход

if (al <= e)or(z_<0) then break;

while (true) do

{3} begin

//отражение относительно 0

kx:= 2/3*(x2+x1+x3)-x0; ky:= 2/3*(y2+y1+y3)-y0; kz:= 2/3*(z2+z1+z3)-z0;

ZZ(x0, y0,z0); zz1:=z_; ZZ(kx, ky,kz); zz2:=z_;

if (z_<0) then goto l;

if (zz1 < zz2)and(kx>0)and(ky>0)and(kz>0)and((6*kx*ky+4*kz*kx)<=s) then begin x0:= kx; y0:= ky; z0:= kz; inc(iter); break; end;

//отражение относительно 1

kx:= 2/3*(x2+x0+x3)-x1; ky:= 2/3*(y2+y0+y3)-y1; kz:= 2/3*(z2+z0+z3)-z1;

ZZ(x1, y1,z1); zz1:=z_; ZZ(kx, ky,kz); zz2:=z_; if (z_<0) then goto l;

if (zz1 < zz2)and(kx>0)and(ky>0)and(kz>0)and((6*kx*ky+4*kz*kx)<=s)then begin x1:= kx; y1:= ky; z1:= kz;inc(iter); break; end;

//отражение относительно 2

kx:= 2/3*(x0+x1+x3)-x2; ky:= 2/3*(y0+y1+y3)-y2; kz:= 2/3*(z0+z1+z3)-z2;

ZZ(x2, y2,z2); zz1:=z_; ZZ(kx, ky,kz); zz2:=z_;if (z_<0) then goto l;

if (zz1 < zz2)and(kx>0)and(ky>0)and(kz>0)and((6*kx*ky+4*kz*kx)<=s)then begin x2:= kx; y2:= ky; z2:= kz;inc(iter); break; end;

//отражение относительно 3

kx:= 2/3*(x2+x1+x0)-x3; ky:= 2/3*(y2+y1+y0)-y3; kz:= 2/3*(z2+z1+z0)-z3;

ZZ(x3, y3,z3); zz1:=z_; ZZ(kx, ky,kz); zz2:=z_;if (z_<0) then goto l;

if (zz1 < zz2)and(kx>0)and(ky>0)and(kz>0)and((6*kx*ky+4*kz*kx)<=s) then begin x3:= kx; y3:= ky; z3:= kz;inc(iter); break; end;

//уменьшаем треугольник

al:= al / 2;

x1:= x0 + al; y1:= y0;z1:=z0;

x2:= x0 + al / 2; y2:= y0 + al * Sin(60);z2:=z0;

x3:= x0 + al/ 2; y3:=y0 ;z3:=z0 + al * Sin(60);break; inc(iter);

{3} end;

{2} end;

l: x:= x3; y:= y3;

{1} end;

procedure TForm1.ZZ(x,y,z:real);

begin z_:=(4/3*x*x*z+7/3*y*x*x)-gm*sqr(s-(6*x*y+4*z*x)); end;

procedure TForm1.Button1Click(Sender: TObject);

begin

RichEdit1.Text :='';iter:=0;

gm:=StrToFloat(Edit1.Text);s:=StrToFloat(Edit2.Text);al:=0.05;e:=0.01;

x:=StrToFloat(Edit3.Text);y:=StrToFloat(Edit4.Text);z:=StrToFloat(Edit5.Text);

Sym(x,y,z);ZZ(x3,y3,z3); //считаем значение функции в найденной точке

RichEdit1.Text:= RichEdit1.Text + 'Число итераций '+ IntToStr(iter) + #13;

RichEdit1.Text:= RichEdit1.Text + 'x3= '+ FloatToStr(x1) +'; y3= '+ FloatToStr(y1)+'z3= '+ FloatToStr(z1) + #13;

ZZ(x3,y3,z3); RichEdit1.Text:= RichEdit1.Text +'f(x,y,z)= '+ FloatToStr(z_) +'('+FloatToStr(6*x3*y3+4*z3*x)+')'+ #13;end;end.