Есть чему поучиться!
Экономико-математическое моделирование
ЛАБОРАТОРНАЯ РАБОТА №1
Аналитическое определение экстремума функции одной и нескольких переменных
.1 Функции одной переменной
а) Постановка задачи: исследовать на экстремум функцию вида
y=f(x)=a0+a1*x+a2*x2 (1.1)
б)Необходимое условие экстремума: первая производная функции должна принимать значение, равное нулю.
dQ (x)/ dx=0 (1.2)
Для функции (1.1) получим
d Q (x)/ dx= a1+2a2*x (1.3)*=-a1/2*a2 (1.4)
в)Достаточное условие экстремума : порядок впервые не обращающейся в ноль производной при выполнении условия (1.2) должен быть чётным. Причём, если значение этой производной положительное, то функция Q (x) имеет минимум, а если значение этой производной отрицательное, то функция Q (x) имеет максимум.
Для функции (1.1) имеем:
d Q (x)/ dx =0
d2 Q (x)/ dx2 =2*a2
При выполнении необходимого условия экстремума (1.2) получили порядок впервые не обращающийся в ноль производной четный (он равен 2). Следовательно, функция Q (x) имеет экстремум. Знак второй производной в данном примере определяется знаком коэффициента a2.
1.2Функция многих переменных
а) Постановка задачи: исследовать на экстремум функцию вида:
Q(x1,x2)=a0+a1*x1+a2*x12+a3*x2+a4*x22+ a5* x1*x2. (1.5)
б) Необходимое условие экстремума: частные производные функции должны принимать значения, равные нулю
Q (x1,…, xn)/xi=0, i=1,…,n (1.6)
Для функции (1.5) получим
экстремум функция оптимизация инвестиция
¶ Q (x1, х2)/ ¶x1=0,
2*a2*x1+a5*x2=-a1. (1.7)
¶ Q (x1, х2)/ ¶x2=0,
a5*x+2a4*x2=-a3 (1.8)
По формуле Крамера :
X 1=D1/D, (1.9)
X 2=D2/D (1.10)
D=4*a2*a4-a52 (1.11)
D1=-2*a1*a4+a3*a5 (1.12)
D2=-2*a2*a3+a1*a5 (1.13)
Систему уравнений (1.7) и (1.8) можно решать методом подстановки.
в)Достаточное условие экстремума.
Разложим функцию Q (x1…,xn) в ряд Тейлора относительно (x1*,…,xn*) и с учётом (1.6) получим
DQ(x1*,…,xn*)=0,5ai,j*zi*zj, (1.14)
где=Dxi=xi - |xi*|=Dxj=xj-|xj*|
Выражение (1.14 ) называется квадратичной формой.
Составляется матрица из вторых частных производныхj = ¶Q(x1,…,xn)/ ¶xi¶xj
в виде
(1.15)
Квадратичная форма (1.14) является положительно-определённой, если все диагональные миноры матрицы (1.15) строго положительные. Квадратичная форма (1.14) является отрицательно-определённой, если все нечётные диагональные миноры матрицы (1.15) строго отрицательные, а чётные диагональные миноры матрицы (1.15) строго положительные.
Определение:
Функция Q (x1,…,xn) имеет минимум, если квадратичная форма (1.14) строго положительная.
Функция Q (x1,…,xn) имеет максимум, если квадратичная форма (1.14) строго отрицательная.
Для функции (1.5) имеем:= ¶Q(x1,x2)2/¶x1 2=2*a2
a22= ¶Q(x1,x2)2/¶x2 2=2*a4= ¶Q(x1,x2)2 /¶x1¶x2 = a21= ¶Q(x1,x2)2 /¶x2¶x1 =a5
Ниже приведены значения коэффициентов в уравнениях (1.1)
( таблица 1) и (5) (таблица2)
Таблица 1
Значения коэффициентов уравнения 1№ п/п№ п/п12-20202920-381124-22193018-401236-24183116-421348-26173214-4414510-28163312-4615612-30153410-4816714-3214358-5017816-3413366-5218918-3612374-54191020-3811382-56201122-40103920-20201224-4294022-22191326-4484124-24181428-4674226-26171530-4864328-28161632-5054430-301517 34-5244532-32141836-54346 34-34131938-5624736-36122038-2024838-38112136-2234940-40102234-2445042-4282332-2655144-4462430-2865246-4642528-3075348-4822626-3285450-5022724-34915552-5022822-36105648-502
Таблица 2
№ п/п1100382515120210137352321931023645322184103355541617510434655121661053375610157106328578148107319586139108301059412101092911603111111028126131012111271362291311226146328141132515642715114241665161611523156615171162214671418117211368131911820126912201191911701321120181071142212117972152312216873262412315774272512414675382612513576592712612477810281271131112112912810212201230129911350133113081714314321317181541533132619165163413352017517351344211851836135322195193713622320620381371242162139138125252240138226262341138327372442138428472443138529582444139630682545140732782546141834782547141835782648142836782649143837782750144838782951145839783052146840783153147841793154148842793255149842710335615084371134
Каждому студенту задаётся вариант индивидуального задания для исследования функций одной и нескольких переменных на экстремум соответственно из таблиц 1 и 2.
.3 Пример выполнения задания
Постановка задачи
Исследовать на экстремум функции:
- y=10-48x+16x2;
- y=133+5x1+20x12+17x2+5x22+17x1x2.
Задание 1. Необходимое условие экстремума функции одной переменной:
; x=x*.
y=10-48x+16x2;
Находим производную, приравниваем её к нулю:
;
Достаточное условие экстремума функции одной переменной:
;
.
Функция y=10-48x+16x2 имеет минимум при x=1,5, так как порядок производной, впервые не обращающейся в нуль, четный (второй), и значение этой производной больше нуля.
Ответ: y*=f= -26- минимум функции y=10-48x+16x2.
Задание 2. Необходимое условие экстремума функции многих переменных:
y=f(x1,…,xn).
;=133+5x1+20x12+17x2+5x22+17x1x2;
Выражаем x1 из первого уравнения и подставляем во второе.
Получаем уравнение:
-17* (5+17 x2)/40 + 10 x2 = -17;
,775 x2= -14,875;
;
Подставляя значение x2 в первое уравнение, получаем значение для x1:
;
;
;
Достаточное условие экстремума функций многих переменных
Составляем матрицу из коэффициентов квадратичной формы. Для этого вычисляем вторые частные производные исходной функции:
;
;
;
.
Диагональные миноры матрицы квадратов А:
;
;
Значит, квадратичная форма является положительно определенной, т. е. функция
y=133+5x1+20x12+17x2+5x22+17x1x2
при будет иметь минимум.
ЛАБОРАТОРНАЯ РАБОТА №2
Определение оптимальной долговечности изделия аналитическим методом
.1 Постановка задачи
Решается задача определения оптимальной долговечности (срока службы) изделия, участвующего в технологическом процессе, при котором затраты на планируемом временном отрезке Т при использовании изделия с долговечностью Т¶ будут минимальны.
Суммарные затраты:
(2.1)
(2.2) - стоимость НИОКР ;
(2.3) - стоимость серийного производства 1 изделия с долговечностью Т¶;
(2.4)
стоимость годовой эксплуатации 1 изделия с долговечностью Т¶;
Т¶0 - выбранное базовое значение долговечности;
S1, S2, S3 - стоимости НИОКР 1 изделия, серийного производства 1 изделия, стоимости годовой эксплуатации 1 изделия при долговечности
Т¶ = Т¶0. К1П, К2П, КНО, К1ГЭ - коэффициенты, заданные числа.
2.2 Метод решения задачи
2.2.1 Решение задачи аналитическим методом по заданной целевой функции
Задаются значения следующих коэффициентов:
(2.5)
Подставив выражения (2.2), (2.3), (2.4) в выражение (1) с заданными коэффициентами и значениями T и m, а функция (1) примет вид:
(2.6)
Необходимые условия экстремума:
(2.7)
(2.8)
(2.9)
(2.10)
2.2.2 Решение задачи при аппроксимации исходной целевой функции в ряд Тейлора
Разложим функцию (2.7) в ряд Тейлора, ограничившись вторым порядком производной:
(2.11)
Где .
Где Тд0 - базовое значение Тд. Введем обозначения:
(2.12)
; (2.13)
Подставив (2.12) и (2.13)в (2.11), получим выражение функции одной переменной, приведённое в лабораторной работе №1.
Требуется исследовать её на экстремум и полученное оптимальное значение сравнить со значением, вычисленным по исходной целевой функции, вычислив относительную величину отклонения этих значений.
2.3 Пример выполнения задания
2.3.1 Постановка задачи
Решается задача определения оптимальной долговечности (срока службы) изделия, участвующего в технологическом процессе, при котором затраты на планируемом временном отрезке Т при использовании изделия с долговечностью Т¶ будут минимальны.
Суммарные затраты:
(1)
(2) стоимость НИОКР ;
(3)
стоимость серийного производства 1 изделия с долговечностью Т¶;
(4)
стоимость годовой эксплуатации 1 изделия с долговечностью Т¶;
Т¶0 - выбранное базовое значение долговечности;
S1, S2, S3 - стоимости НИОКР 1 изделия, серийного производства 1 изделия, стоимости годовой эксплуатации 1 изделия при долговечности
Т¶ = Т¶0. К1П, К2П, КНО, К1ГЭ - коэффициенты, заданные числа.
.3.2 Метод решения задачи
Решение задачи аналитическим методом по заданной целевой функции
Заданы следующие значения коэффициентов:
Подставим выражения (2), (3), (4) в выражение (1) с заданными коэффициентами и значениями T и m. Тогда функция (1) примет вид:
(6)
Подставив заданные значения коэффициентов (5) в выражение (6), получим:
(7)
Необходимые условия экстремума:
(8)
Тогда оптимальное значение долговечности (8) будет равно:
(9)
Решение задачи при аппроксимации исходной целевой функции в ряд Тейлора
Разложим функцию (7) в ряд Тейлора, ограничившись вторым порядком производной:
(10)
Где (11)
Из выражения (10) с учетом (11) получим:
(12)
- Необходимое условие экстремума:
- Допустимое условие экстремума:
Так как порядок, впервые не обращающейся в нуль производной равен 2 (четный), то функция имеет экстремум, а так как эта производная положительная, то функция имеет минимум.
2.3.3 Результаты расчетов
Решая аналитическим методом данную задачу, мы получили оптимальное значение долговечности, равное . Оно отличается от оптимального значения долговечности, найденного методом разложения в ряд Тейлора, равного на 25%.
ЛАБОРАТОРНАЯ РАБОТА №3
Решение одномерной задачи статической оптимизации численными методами
.1 Постановка задачи
Требуется найти экстремум функции , где .
Левая граница интервала поиска: .
Правая граница интервала поиска: ,
где - точка экстремума функции , найденная аналитическим методом.
Погрешность нахождения экстремума
При этом нахождение экстремума должно реализовываться следующими методами:
половинного деления;
"золотого" сечения;
с использованием чисел Фибоначчи.
.2 Метод половинного деления.
Определяется интервал поиска (a,b).
Задается точность ? - точность местоположения экстремума от переменной u.
Сущность метода (на примере нахождения минимума):
Исходный интервал делится пополам, т.е. (b1 + a1)/2 = u1
Далее делаются шаги вправо и влево на величину ? = (0.1, … ,0.5) ?.
- ?, находим значения Q()
+ ?, находим значения Q()
Затем сравниваем эти значения.
Если , рассматриваем новый отрезок (ai+1 = ; bi+1 = bi).
Затем сравниваем эти значения.
Далее процедура сокращения интервала поиска повторяется.
Процедуры сокращения интервала поиска заканчиваются когда (b - a)·(1/2)i = < ? .
Число обращений функций - N = i·2.
3.3 Метод золотого сечения
В этом методе отрезок [a,b] разбивается точками x1 и x2, используя пропорцию золотого сечения:
|ab|/ |x1b| = |x1b| / |ax1| (для точки x1)
|ab|/ |ax2| = |ax2| / |x2b| (для точки x2)
К золотому сечению приводит деление отрезка на две неравные части таким образом, что отношение длины всего отрезка к длине его большей части равно отношению длины большей части к длине меньшей части.
Здесь и далее x1 будет левее x2, a |ax1| = |x2b|,
т.е. если известна x1, то x2 = b - (x1 - a),
а если известна x2, то x1 = a + (b - x2).
Находим значение функции в этих точках. В случае минимизации функции при условии F1 F2 отрезок [a,b] сужается до длины с новыми границами: a1 = a, b1 = x2, а при условии F2 < F1 отрезок [a,b] сужается до длины с новыми границами: a1 = х1, b1 = b.
В случае максимизации функции при условии F1 F2 отрезок [a,b] сужается до длины с новыми границами: a1 = a, b1 = x2, а при условии F2 F1 отрезок [a,b] сужается до длины с новыми границами: a1 = х1, b1 = b.
Далее процедура использования стратегии золотого сечения повторяется для большего отрезка в полученном интервале поиска. Поиск заканчивается, когда bi - ai <= ?, где i - число сокращений интервала поиска.
Число вычислений функции N = i+1
3.4 Метод с использованием чисел Фибоначчи.
Используются следующие числа Фибоначчи, формируемые по алгоритму: F0 = F1 = 1; для Fi = Fi-2 + Fi-1 для i=2,3….. Например
I0,1,2,3,4,5,6,7…..
Fi1,1,2,3,5,8,13,2…..
Вначале вычисляют число N=(b-a)/. Затем в ряду чисел Фибоначчи находят пару соседних чисел, удовлетворяющих условию. Далее вычисляют X1=а + dx*F(s-2) и X2 = а + dx*F(s-2-1) и значение функции в этих точках: Q(X1) и Q(X2). Далее движение в сторону экстремума выполняется по общей формуле: Xi+1 = Xt ± sign (Xi-Xi-1)·sign |Q(Xi)-Q(Xi-1)| ?·Fs-2-i , i = 2,3,… Здесь знак +ставится в задаче максимизации функции, а знак - ставится в задаче минимизации функции.
Xt есть значение X, при котором получено наилучшее значение функции на предыдущих шагах. Sign (z) есть знаковая функция. Sign (z)=-1, если z0, Sign (z)=1, если z >0.
Признаком окончания поиска является условие, когда весь ряд чисел Фибоначчи будет исчерпан.Число вычислений функции при этом равно s-1.
3.5 Пример выполнения работы
Постановка задачи.
Требуется найти экстремум функции , где
Левая граница интервала поиска: .
Правая граница интервала поиска: , так как - точка экстремума функции , найденная аналитическим методом.
Погрешность нахождения экстремума
.5.1 Поиск экстремума функции методом половинного деления
График функции
Левая граница отрезка [A, B]: А = 1,05
Правая граница отрезка [А, B]: B = 2,1
Требуемая точность: E = 0,01
б = 0.05*E = 0,0005
Шаг N: 1
Текущий отрезок [A, B]:= 1,05= 2,1A = 1,05
Так как отрезок B-A > E, то отрезок сужается.
Отложим по обе стороны от середины текущего отрезка [A, B] точки= (A + B)/2 - б = 1,5745= (A + B)/2 + б = 1,5755
и отстоящие от середины отрезка [A, B] на расстояние б
Вычислим значения функции в точках U1 и U2= F(U1) = -25,911196= F(U2) = -25,908796
Так как Q1 <= Q2, то далее сужается отрезок [A, B]:
A = A = 1,05= U2 = 1,5755
Шаг N: 2
Текущий отрезок [A, B]:
A = 1,05= 1,5755A = 0,5255
Так как отрезок B-A > E, то отрезок сужается.
Отложим по обе стороны от середины текущего отрезка [A, B] точки= (A + B)/2 - б = 1,31225= (A + B)/2 + б = 1,31325
и отстоящие от середины отрезка [A, B] на расстояние б
Вычислим значения функции в точках U1 и U2= F(U1) = -25,435999= F(U2) = -25,441991
Так как Q1 > Q2, то далее сужается отрезок [A, B]:= U1 = 1,31225= B = 1,5755
……………………..
Шаг N: 7
Текущий отрезок [A, B]:= 1,492546875= 1,5099375A = 0,017390625
Так как отрезок B-A > E, то отрезок сужается.
Отложим по обе стороны от середины текущего отрезка [A, B] точки= (A + B)/2 - б = 1,5007421875= (A + B)/2 + б = 1,5017421875
и отстоящие от середины отрезка [A, B] на расстояние б
Вычислим значения функции в точках U1 и U2= F(U1) = -25,9999911865234= F(U2) = -25,9999514365234
Так как Q1 <= Q2, то далее сужается отрезок [A, B]:
A = A = 1,492546875= U2 = 1,5017421875
Шаг N: 8
Текущий отрезок [A, B]:
A = 1,492546875= 1,5017421875A = 0,00919531249999994
Так как отрезок B-A < E, то требуемая точность достигнута.
Точка минимума: Xmin = B = 1,5017421875
Значение функции в точке Xmin: F(Xmin) = -25,9999514365234
Поиск окончен.
.5.2 Поиск экстремума функции методом Золотого сечения
Описание алгоритма поиска
Исследуемая фунция: 10-48*x+16*x*x
Левая граница отрезка [A, B]: А = 1,05
Правая граница отрезка [А, B]: B = 2,1
Требуемая точность: E = 0,01
Разобьем текущий отрезок [A; B] точками X1 и X2, используя пропорцию золотого сечения:
|AB| / |X1 B| = |X1 B| / |A X1| (для точки X1)
|AB| / |A X2| = |A X2| / |X2 B| (для точки X2)
Здесь и далее X1 будет левее X2, a |AX1| = |X2B|,
т.е. если известна X1, то X2 = B - (X1 - A),
а если известна X2, то X1 = A + (B - X2)
Шаг N: 1
Текущий отрезок [A, B]:= 1,05= 2,1A = 1,05
Так как отрезок B-A > E, то отрезок сужается.= 1,449, F1 = F(X1) = -25,958384= 1,701, F2 = F(X2) = -25,353584
Так как F1 <= F2, то далее сужается отрезок [A, B]:
A = A = 1,05= X2 = 1,701
Шаг N: 2
Текущий отрезок [A, B]:
A = 1,05= 1,701A = 0,651
Так как отрезок B-A > E, то отрезок сужается.= 1,29738, F1 = F(X1) = -25,3431221696= 1,45362, F2 = F(X2) = -25,9655823296
Так как F1 > F2, то далее сужается отрезок [A, B]:= X1 = 1,29738= B = 1,701
………………..
Шаг N: 10
Текущий отрезок [A, B]:= 1,49602107020964= 1,51023501108321A = 0,0142139408735766
Так как отрезок B-A > E, то отрезок сужается.= 1,5014223677416, F1 = F(X1) = -25,9999676299201= 1,50483371355125, F2 = F(X2) = -25,9996261634129
Так как F1 <= F2, то далее сужается отрезок [A, B]:
A = A = 1,49602107020964= X2 = 1,50483371355125
Шаг N: 11
Текущий отрезок [A, B]:
A = 1,49602107020964= 1,50483371355125A = 0,00881264334161758
Так как отрезок B-A < E, то требуемая точность достигнута.
Точка минимума: Xmin = A = 1,49602107020964
Значение функции в точке Xmin: F(Xmin) = -25,9997466898836
Поиск окончен.
.5.3 Поиск экстремума функции методом с использованием чисел Фибоначчи
Исследуемая функция: 10-48*x+16*x*x
Левая граница отрезка [A, B]: А = 1,05
Правая граница отрезка [А, B]: B = 2,1
Требуемая точность: E = 0,01= (B-A)/E = 105
F(s-1) = 89 < (N = 105) <= F(s) = 144
Элементарный шаг - delta (dm) = (B-A)/F(s) = 0,00729166666666667
Подготовительный шаг
Вычислим значение функции на левой границе отрезка.= A = 1,05= Q(X1) = -22,76
Будут использоваться следующие числа Фибоначчи
Шаг N: 1
Шаг делается из точки X1 = 1,05,
при этом Q(X1) = -22,76 вычислено шагом раньше.
Величина шага dx = (dm = 0,00729166666666667) * (F(s-2-0)=55) = 0,401041666666667= X1 + dx = 1,05 + 0,401041666666667=1,45104166666667(X2) = -25,9616493055556
Шаг оказался удачным, так как (Q2 = -25,9616493055556) <= (R1 = -22,76)
Точка X2 оказалась удачнее точки X1 из которой делали шаг.
Следующий шаг будет сделан из удачной точки полученной на этом шаге:= X2 = 1,45104166666667,
При этом будет использовано Q1 = Q(X1=X2) = -25,9616493055556 уже вычисленное на проделанном шаге
Шаг N: 2
Шаг делается из точки X1 = 1,45104166666667,
при этом Q(X1) = -25,9616493055556 вычислено шагом раньше.
Величина шага dx = (dm = 0,00729166666666667) * (F(s-2-1)=34) = 0,247916666666667= X1 + dx = 1,45104166666667 + 0,247916666666667=1,69895833333333(X2) = -25,3666493055556
Шаг оказался неудачным, так как (Q2 = -25,3666493055556) > (R1 = -25,9616493055556)
Точка X2 оказалась хуже точки X1 из которой делали шаг.
Следующий шаг будет сделан из более удачной точки X1,
но в противоположном направлении, т.е. знак приращения будет противоположен
знаку приращения на неудачном шаге.
При этом будет использовано Q1 = Q(X1) = -25,9616493055556 уже вычисленное шагом ранее
Шаг N: 8
Шаг делается из точки X1 = 1,509375,
при этом Q(X1) = -25,99859375 вычислено шагом раньше.
Величина шага dx = (dm = 0,00729166666666667) * (F(s-2-7)=2) = 0,0145833333333333= X1 + dx = 1,509375 + 0,0145833333333333=1,52395833333333(X2) = -25,9908159722222
Шаг оказался неудачным, так как (Q2 = -25,9908159722222) > (R1 = -25,99859375)
Точка X2 оказалась хуже точки X1 из которой делали шаг.
Следующий шаг будет сделан из более удачной точки X1,
но в противоположном направлении, т.е. знак приращения будет противоположен
знаку приращения на неудачном шаге.
При этом будет использовано Q1 = Q(X1) = -25,99859375 уже вычисленное шагом ранее
Шаг N: 9
Шаг делается из точки X1 = 1,509375,
при этом Q(X1) = -25,99859375 вычислено шагом раньше.
Величина шага dx = (dm = 0,00729166666666667) * (F(s-2-8)=1) = 0,00729166666666667= X1 - dx = 1,509375 - 0,00729166666666667=1,50208333333333(X2) = -25,9999305555556
Шаг оказался удачным, так как (Q2 = -25,9999305555556) <= (R1 = -25,99859375)
Точка X2 оказалась удачнее точки X1 из которой делали шаг.
Следующий шаг будет сделан из удачной точки полученной на этом шаге:= X2 = 1,50208333333333,
При этом будет использовано Q1 = Q(X1=X2) = -25,9999305555556 уже вычисленное на проделанном шаге.
Все числа Фибоначчи исчерпаны. Итерации завершены.= -25,9999305555556 Q2 = -25,9999305555556
Так как Q2 <= Q1 то предыдущий шаг оказался удачным
В качестве точки минимума принимается точка X2, которая оказалась лучше точки X1= 1,50208333333333.
Сводная таблица результатов
№ методаНаименование методаX*Y(X*)e Число обращений к функции1Аналитический метод1,5-26--2Метод половинного деления1,50-25,990,01143Метод золотого сечения1,49-25,990,01114Метод с использованием чисел Фибоначчи1,50-25,990,0110
3.6 Вывод
Сравнивая результаты, полученные при вычислении экстремума функции численными методами с результатами аналитического метода, видим, что численные методы дают приблизительное значение экстремума функции, близкое к точному. Наиболее эффективным методом поиска оказался метод с использованием чисел Фибоначчи
ЛАБОРАТОРНАЯ РАБОТА №4
Решение многомерной задачи статической оптимизации численными методами
.. Постановка задачи:
Задана функция:
(4.1)
конкретным выражением
,
где
Границы интервалов поиска:
Предельные значения производных:
,
Начальные шаги по координатным осям:
ДельтаX
ДельтаY
Базовый шаг в направлении градиента: h0
Заданное предельное значение модуля градиента: E
Найти экстремум функции методами:
Гаусса-Зейделя, градиента, наискорейшего спуска.
Требуется найти минимум функции.
4.2 Метод Гаусса-Зейделя
Описание алгоритма поиска
В этом методе последовательность переменных определяется рядом
Алгоритм метода
Выбирается исходная точка поиска
Вычисляется новая точка , где - величина рабочего шага.
Знак "+" ставится а знак в задаче максимизации,
а знак " - " ставится в задаче минимизации.
Вычисляется величина
Далее вычисляем
Для произвольной переменной
Останов поиска в направлении выбранной переменной производится , если функция перестает улучшаться , то есть при , . Затем изменяют следующую по порядку переменную до выполнения аналогичного условия.
Поиск заканчивается при выполнении условия:
,
4.3 Метод градиента
Алгоритм метода
Постановка задачи:
Формула градиента функции:
где - направляющие вектора по оси ., где - точка приложения
Признак окончания поиска:
Направление определяется направляющими косинусами:
где .
Алгоритм поиска методом градиента:
1.Выбирается исходная точка поиска в области допустимых управлений (решений):
.
В ней вычисляется значение функции, направляющие косинусы и оценки частных производных .
Перейдем от частных производных к отношению приращений:
. Направляющие косинусы вычисляются по следующей формуле:
где .
1.Формируется шаг в направлении градиента следующим образом:
, где
Знак "+" ставится в задачах максимизации,
"-" - в задачах минимизации.
Существуют различные способы формирования шага. Широко распространен способ формирования шага пропорционально модулю градиента.
.
- начальная величина шага.
2.Алгоритм движения вдоль j-й координатной оси:
При каждом шаге в направлении градиента вычисляется целевая функция в полученной точке и частные производные по каждой из переменных .
3.Проверяется выполнение признака окончания поиска:
, где малая величина.
4.4 Методом наискорейшего спуска
Алгоритма метода
Метод является синтезированным из двух описанных выше. В данном методе движение из исходной точки происходит в направлении градиента n до выполнения условия , . Из точки, в которой выполнено это условие, движение осуществляется в направлении градиента, т.е. перпендикулярно предыдущему направлению. Поиск заканчивается по тому же условию, что и в двух предыдущих.
4.5 Пример выполнения задания
Исследуемая функция F(X, Y) = 100+38*x+2*x*x+51*y+51*y*y+20*x*y
Интервал по X = [-374,86848,-339,16672]
Интервал по Y = [66,0288,72,9792]
Требуемая точность: E = 0,01
Начальный шаг h0 = 0,01
Сигма sigX = 0,0001
Сигма sigY = 0,0001
ДельтаX = 1E-5
ДельтаY = 1E-5
Координаты стартовой точки:
X начальное = -365,94304
Y начальное = 67,7664
.5.1 Метод Гаусса-Зейделя
Значение функции в текущей точке Fn(x, y) F0 = -4287,44717127684
Значение частной производной по x = -70,4441335983574
Значение частной производной по y = -355,687492992729
< . Начинаем поиск вдоль по оси оy
Вычислим новую точку y = y - sign()*h = 67,7689
Шаг 1
Значение функции в текущей точке F(x, y) = -4288,33607252681
Значение частной производной по y = -355,432491051033
Значение функции лучше чем предыдущее. Fn = -4288,33607252681 F0 = -4287,44717127684
Вычислим новую точку y = y - sign() * h = 67,7714
F0 = Fn = -4288,33607252681
Шаг 2
Значение функции в текущей точке F(x, y) = -4289,22433627688
Значение частной производной по y = -355,177489109337
Значение функции лучше чем предыдущее. Fn = -4289,22433627688 F0 = -4288,33607252681
Вычислим новую точку y = y - sign() * h = 67,7739
F0 = Fn = -4289,22433627688
Шаг 3
Значение функции в текущей точке F(x, y) = -4290,1119625268
Значение частной производной по y = -354,922498809174
Значение функции лучше чем предыдущее. Fn = -4290,1119625268 F0 = -4289,22433627688
Вычислим новую точку y = y - sign() * h = 67,7764
F0 = Fn = -4290,1119625268
Шаг 2970
Значение функции в текущей точке Fn(x, y) = -4910,409048045
Значение частной производной по x = -0,0100815668702126
Значение функции Fn -4910,409048045 лучше чем предыдущее F0 = -4910,40903229825
Вычислим новую точку x = x - sign() * h = -359,920774375
F0 = Fn = -4910,409048045
Шаг 2971
Значение функции в текущей точке Fn(x, y) = -4910,40905754192
Значение частной производной по x = -0,00507570803165436
Значение функции Fn -4910,40905754192 лучше чем предыдущее F0 = -4910,409048045
Вычислим новую точку x = x - sign() * h = -359,919524375
F0 = Fn = -4910,40905754192
Шаг 2972
Значение функции в текущей точке Fn(x, y) = -4910,40906078881
Значение частной производной по x = -8,14907252788544E-5
Значение функции Fn -4910,40906078881 лучше чем предыдущее F0 = -4910,40905754192
Частная производная = -8,14907252788544E-5 меньше sigx = 0,0001
Проверка на достижения условия выхода
Корень суммы квадратов частных производных A = 0,903524739755089
Так как A (0,903524739755089) <= E (0,1) то поиск окончен.
Экстремум найден.
Xmin = -359,919524375 Ymin = 70,0838999999962
Значение функции в точке экстремума F(Xmin, Ymin) = -4910,40906078881
Экстремум найден.
Xmin = -357,148864228791 Ymin = 69,5297901322648
Значение функции в точке экстремума F(Xmin, Ymin) = -4910,74911253306
График функции
4.5.2 Метод градиента
Шаг 1
Делаем шаг из точки X = -365,94304 Y = 67,7664
Функция F(X, Y) в данной точке = -4287,44717127684
Частная производная (по X) = -70,4439589753747
Частная производная (по Y) = -355,682899826206
Проверка на выход:
Корень суммы квадратов частных производных A = 362,591611299687
Так как A (362,591611299687) > E (0,1) то продолжаем поиск.
Вычислим новые точки X и Y
Xновое = X - h0*= -365,238600410246
Yновое = Y - h0*= 71,3232289982621
Текущая точка по X приравнивается к новой вычисленной точке по X = -365,238600410246
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 71,3232289982621
Шаг 2
Делаем шаг из точки X = -365,238600410246 Y = 71,3232289982621
Функция F(X, Y) в данной точке = -4905,88566703646
Частная производная (по X) = 3,51038004737347
Частная производная (по Y) = 21,2024501524866
Проверка на выход:
Корень суммы квадратов частных производных A = 21,4910832799483
Так как A (21,4910832799483) > E (0,1) то продолжаем поиск.
Вычислим новые точки X и Y
Xновое = X - h0*= -365,27370421072
Yновое = Y - h0*= 71,1112044967372
Текущая точка по X приравнивается к новой вычисленной точке по X = -365,27370421072
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 71,1112044967372
Шаг 3
Делаем шаг из точки X = -365,27370421072 Y = 71,1112044967372
Функция F(X, Y) в данной точке = -4908,05924940115
Частная производная (по X) = -0,870526419021189
Частная производная (по Y) = -1,12612557131797
Проверка на выход:
Корень суммы квадратов частных производных A = 1,42336750299776
Так как A (1,42336750299776) > E (0,1) то продолжаем поиск.
Вычислим новые точки X и Y
X новое = X - h0*= -365,26499894653
Y новое = Y - h0*= 71,1224657524504
Текущая точка по X приравнивается к новой вычисленной точке по X = -365,26499894653
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 71,1224657524504
Шаг 2453
Делаем шаг из точки X = -358,290583220967 Y = 69,7531945287842
Функция F(X, Y) в данной точке = -4910,68468116963
Частная производная (по X) = -0,098241725936532
Частная производная (по Y) = 0,0192777952179313
Проверка на выход:
Корень суммы квадратов частных производных A = 0,100115284065187
Так как A (0,100115284065187) > E (0,1) то продолжаем поиск.
Вычислим новые точки X и Y
X новое = X - h0*= -358,289600803708
Y новое = Y - h0*= 69,7530017508321
Текущая точка по X приравнивается к новой вычисленной точке по X = -358,289600803708
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 69,7530017508321
Шаг 2454
Делаем шаг из точки X = -358,289600803708 Y = 69,7530017508321
Функция F(X, Y) в данной точке = -4910,68478057638
Частная производная (по X) = -0,0981672201305628
Частная производная (по Y) = 0,0192632433027029
Проверка на выход:
Корень суммы квадратов частных производных A = 0,100039370503325
Так как A (0,100039370503325) > E (0,1) то продолжаем поиск.
Вычислим новые точки X и Y
X новое = X - h0*= -358,288619131506
Y новое = Y - h0*= 69,752809118399
Текущая точка по X приравнивается к новой вычисленной точке по X = -358,288619131506
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 69,752809118399
Шаг 2455
Делаем шаг из точки X = -358,288619131506 Y = 69,752809118399
Функция F(X, Y) в данной точке = -4910,68487983197
Частная производная (по X) = -0,0980938784778118
Частная производная (по Y) = 0,0192475272342563
Проверка на выход:
Корень суммы квадратов частных производных A = 0,0999643751516166
Так как A (0,0999643751516166) <= E (0,1) то поиск окончен.
Экстремум найден.
Xmin = -358,288619131506 Ymin = 69,752809118399
Значение функции в точке экстремума F(Xmin, Ymin) = -4910,68487983197
График функции
4.5.3 Метод наискорейшего спуска
Шаг 1
Делаем шаг из точки X = -365,94304 Y = 67,7664
Функция F(X, Y) в данной точке = -4287,44717127684
Текущий шаг h = 0,1
Значение функции в текущей точке F(X, Y) = -4287,44717127684
Значение частной производной по x = -70,4439589753747
Значение частной производной по y = -355,682899826206
Проверка на выход.
Корень суммы квадратов частных производных A = 362,591611299687
Так как A (362,591611299687) > E (0,1) то продолжаем поиск.
Меняем направление. Вычислим направляющие косинусы.
Направляющий косинус по X cosX = -0,194279064324938
Направляющий косинус по Y cosY = -0,980946300856997
Значение частной производной по x = -70,4439589753747
Значение частной производной по y = -355,682899826206
Вычислим новую точку поиска.
X новое = X - sign(dfdx)*|cosX|*h = -365,923612093567
Y новое = Y - sign(dfdy)*|cosY|*h = 67,8644946300857
Значение функции в новой точке F(X,Y) = -4323,1772158799
Проверка удачности шага.
Текущая точка по X приравнивается к новой вычисленной точке по X = -365,923612093567
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 67,8644946300857
Шаг 2
Делаем шаг из точки X = -365,923612093567 Y = 67,8644946300857
Функция F(X, Y) в данной точке = -4323,1772158799
Значение частной производной по x = -68,4043555520475
Значение частной производной по y = -345,288689713925
Вычислим новую точку поиска.
X новое = X - sign()*|cosX|*h = -365,904184187135
Y новое = Y - sign()*|cosY|*h = 67,9625892601714
Значение функции в новой точке F(X,Y) = -4357,84801901889
Проверка удачности шага.
Текущая точка по X приравнивается к новой вычисленной точке по X = -365,904184187135
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 67,9625892601714
Шаг 3
Делаем шаг из точки X = -365,904184187135 Y = 67,9625892601714
Функция F(X, Y) в данной точке = -4357,84801901889
Значение частной производной по x = -66,3647521287203
Значение частной производной по y = -334,894479601644
Вычислим новую точку поиска.
X новое = X - sign()*|cosX|*h = -365,884756280703
Y новое = Y - sign()*|cosY|*h = 68,0606838902571
Значение функции в новой точке F(X,Y) = -4391,45958069392
Проверка удачности шага.
Текущая точка по X приравнивается к новой вычисленной точке по X = -365,884756280703
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 68,0606838902571
Шаг 11449
Делаем шаг из точки X = -357,150410737705 Y = 69,5300131040769
Функция F(X, Y) в данной точке = -4910,74909897527
Значение частной производной по x = -0,00118045136332512
Значение частной производной по y = 0,0582222128286958
Вычислим новую точку поиска.
X новое = X - sign()*|cosX|*h = -357,148864228791
Y новое = Y - sign()*|cosY|*h = 69,5297901322648
Значение функции в новой точке F(X,Y) = -4910,74911253306
Проверка удачности шага.
Текущая точка по X приравнивается к новой вычисленной точке по X = -357,148864228791
Текущая точка по Y приравнивается к новой вычисленной точке по Y = 69,5297901322648
Шаг 11450
Делаем шаг из точки X = -357,148864228791 Y = 69,5297901322648
Функция F(X, Y) в данной точке = -4910,74911253306
Значение частной производной по x = 0,00054540578275919
Значение частной производной по y = 0,066409120336175
Вычислим новую точку поиска.
X новое = X - sign()*|cosX|*h = -357,150410737705
Y новое = Y - sign()*|cosY|*h = 69,5295671604527
Значение функции в новой точке F(X,Y) = -4910,74911252246
Проверка удачности шага.
Так как F(Xновое, Yновое) = -4910,74911252246>=F(X, Y) = -4910,74911253306
то делим шаг пополам.
Новый шаг h = 0,00078125
Так как текущее значение h < sig (0,001), то приравниваем h к h0.
Шаг 11451
Делаем шаг из точки X = -357,148864228791 Y = 69,5297901322648
Функция F(X, Y) в данной точке = -4910,74911253306
Текущий шаг h = 0,1
Значение функции в текущей точке F(X, Y) = -4910,74911253306
Значение частной производной по x = 0,00054540578275919
Значение частной производной по y = 0,066409120336175
Проверка на выход.
Корень суммы квадратов частных производных A = 0,0664113599566553
Так как A (0,0664113599566553) <= E (0,1) то поиск окончен.
График функции
Итоговая таблица
№ методаНазвание методаКол-во шаговКол-во вычислений целевой функции1Аналитический метод-357,0169,50-4910,74---2Метод Гаусса - Зейделя-359,9170,08-4910,400,1297233283Метод Градиента-358,2869,75-4910,680,1245573664Метод Наискорейшего спуска-357,1469,52-4910,740,11145145817
Вывод: Сравнив результаты, полученные при вычислении экстремума функции численными методами с результатами аналитического метода получили, что численные методы дают приблизительное значение экстремума функции, близкое к точному. Из численных методов самым эффективным оказался метод Градиента.
ЛАБОРАТОРНАЯ РАБОТА №5
Решение задачи оптимального вложения инвестиций, как задачи распределения, методом динамического программирования
.1 Постановка задачи
Задан объем инвестиций А = 5. Требуется оптимально распределить этот объем на три мероприятия. Эффективность каждого мероприятия от вложенных инвестиций оценивается величиной прибыли qi(Ui).
(5.1)
А - общий объем инвестиций;
Ui = 0, 1, 2, 3, 4, 5 - объем инвестиций в каждое мероприятие;
qi(Ui) - эффективность вложенных средств, которая выражается в виде составляющей эффекта вложенных средств (результат вложения средств в i-тое мероприятие - прибыль).
Нам даны следующие значения:
q1(U1) = 1,2*0,2*U1 =0,24 U1.
U2012345q2 00,090,4050,450,540,63U3012345q3 00,550,7150,8250,880,9355.2 Метод решения
Данную задачу следует решать методом динамического программирования.
В основе динамического программирования лежит принцип оптимальности: Оптимальная стратегия обладает тем свойством, что каковы бы ни были первоначальное состояние и управление, при которых система пришла в это состояние, последующее решение (управления) должны быть оптимальными по отношению к этому состоянию.
Сущность динамического программирования состоит в том, что задача большой размерности разбивается на ряд подзадач меньшей размерности (проводится декомпозиция). Это достигается разбивкой всего процесса на несколько стадий. Процесс решения задачи выполняется пошагово. На первом шаге вычисляют значение составляющей критерия оптимальности от одной стадии. На втором и последующих шагах оптимизируется суммарное приращение критерия оптимальности от вновь вводимой стадии процесса вместе с полученным результатом оптимизации на предыдущем шаге. Математически это представляется следующим образом.
Пусть критерий оптимальности имеет вид
(5.2)
где Si,Ui - соответственно состояние и управление i-ой стадии процесса. (i=1,2,…,N).
Последующие состояния определяются выражениями:
(5.3)
Первый шаг оптимизации опишем выражением
Второй и последующий шаги описываются рекуррентным уравнением
(5.4)
Выражения (5.3) и (5.4) называются функциональными уравнениями динамического программирования.
Пример выполнения задания:
шаг.
шаг.
шаг.
5.3 Результаты расчетов
1 шагU1S1000,240,240,480,480,720,720,960,961,21,22 шагS2=X0000000010,240,240,2410,1000,10020,480,480,4810,0910,240,3320,405000,4050030,720,720,7210,0920,480,5720,40510,240,64530,45000,450040,960,960,9610,0930,720,8120,40520,480,88530,4510,240,6940,54000,540051,21,21,210,0940,961,0520,40530,721,12530,4520,480,9340,5410,240,7850,63000,63
3 шагS3=A0051,21,21,5110,5540,961,5120,71530,721,43530,82520,481,30540,8810,241,1250,935000,935
В результате проделанных расчетов получено оптимальное распределение инвестиций объемом в 5 единиц на три мероприятия следующим образом:
Таким образом, следует в первое мероприятие вложить 4 единицы инвестиций и в третье - 1 единицу.
ЛАБОРАТОРНАЯ РАБОТА №6
Решение классической задачи синтеза оптимального управления с использованием принципа максимума
.1 Постановка задачи
Требуется минимизировать время перехода системы из начального состояния y1(0)=ya0 в конечное y(tk)=yk.
Система описывается дифференциальным уравнением:
(6.1)
где y(t) - фазовая координата, а U(t) - управляющее воздействие
6.2 Постановка задачи в математической форме.
(6.2)
(6.3)
(6.4)
(6.5а)
(6.5б)
(6.5в)
(6.5г)
(6.6)
Область допустимых управлений удовлетворяет условию (6.6).
Требуется для заданных начальных [(6.5а) и (6.5б)] и конечных [(6.5в) и (6.5г)] граничных условий состояния системы:
1.построить траекторию движения системы в фазовом пространстве, обеспечивающую минимальное время перехода из начального в конечное;
2.траекторию неоптимального движения системы в фазовом пространстве;
.построить графики изменения фазовых координат системы во времени при переходе системы из начального состояния в конечное за минимальное время;
.построить графики изменения оптимального управления во времени, обеспечивающего минимальное время перехода из исходного состояния в конечное;
.построить графики изменения фазовых координат системы во времени при переходе системы из исходного состояния системы в конечное;
.построить графики изменения неоптимального управляющего воздействия.
.3 Решение задачи принципом максимума
(6.7)
В данной задаче:
(6.8)
С учетом конкретного математического описания системы (6.3) и (6.4) функция Гамильтона принимает вид:
(6.9)
Система дифференциальных уравнений, сопряженных системе (6.3) и (6.4) имеет вид:
(6.10)
(6.11)
Согласно принципу максимума при оптимальном управлении функция Гамильтона (9) принимает максимальное значение. Следовательно, для поиска оптимального управления Uопт(t), нужно максимизировать функцию Гамильтона:
(6.12)
Так как функция Гамильтона (6.9) линейно зависит от управления U(t), то оптимальное управление будет принимать значения: Uопт(t)=-1, если множитель при U(t) в (6.9) имеет знак «-»; и Uопт(t)=1, если множитель при U(t) в (6.9) имеет знак «+», т.е.: Uопт(t)=sign?2(t), где sign?2(t) - знаковая функция.
В соответствии с алгоритмом решения задачи принципом максимума найденное выражение оптимального управления подставляют в систему сопряженных уравнений (6.3), (6.4), (6.10), (6.11).
Решение этой системы при заданных граничных условиях (6.5а), (6.5б), (6.5в) и (6.5г) можно получить в фазовом пространстве
(6.14)
или во временной области
(6.15)
(6.16)
При решении системы уравнений (6.3), (6.4), (6.10), (6.11) в фазовом пространстве оптимальное управление является функцией фазовых координат. При решении системы уравнений (6.3), (6.4), (6.10), (6.11) во временной области оптимальное управление является функцией времени. Покажем это: рассмотрим систему уравнений (6.10) и (6.11).
Из (6.10) следует:
(6.17)
Подставим (6.17) в (6.11):
(6.18)
и проинтегрируем (6.18)
(6.19)
Заметим, что функция (6.19) линейна. Следовательно, ?2(t) только один раз может сменить знак (рис. 1).
Рис. 1. Графики функций ?2(t) и Uопт(t)
В соответствии с (6.13) оптимальное управление тоже только один раз меняет знак и может переключаться со значения (-1) на (+1) или со значения (+1) на (-1).
.4 Построение фазовых траекторий.
Построим фазовые траектории движения системы для случая Uопт(t)=+1. При этом система уравнений (3) и (4) описывающая рассматриваемую систему, примет вид:
(6.20)
(6.21)
Поделим (6.20) на (6.21):
(6.22)
и проинтегрируем (6.22):
(6.23)
(6.24)
Уравнение (6.24) описывает движение системы в фазовом пространстве при U(t)=1. Графики функции (6.24) для различных y0 приведены на рис.2. Они представляют собой параболы с ветвями симметричными относительно оси y1(t), и вершинами слева.
Рис. 2. Графики движения системы в фазовом пространстве в зависимости от начальных условий и функции управления.
Построим фазовые траектории движения системы для случая Uопт(t)=-1. При этом система уравнений (3) и (4), описывающая рассматриваемую систему, примет вид:
(6.25)
(6.26)
Решение системы уравнений (6.25) и (6.26) найдем аналогично предыдущему случаю:
(6.27)
Графики функций для различных y0 приведены на рис. 2. Они представляют собой параболы с ветвями, симметричными относительно оси y1(t), и вершинами справа.
Заметим, что фазовая координата y2(t) является скоростью изменения фазовой координаты y1(t). Поэтому в положительной полуплоскости относительно оси y1(t) (т.е. для y2(t)>0) координата y1(t) может только возрастать, а в отрицательной полуплоскости относительно оси y1(t) (т.е. для y2(t)<0) координата y1(t) может только убывать.
Геометрически задача оптимального управления сводится к тому, что из заданной исходной точки фазового пространства с координатами (y10, y20) нужно перейти в заданную конечную точку с координатами (y1k, y2k) и при этом осуществить лишь один переход с одной траектории, проходящей через точку с координатами (y10, y20) на траекторию с координатами (y1k, y2k). Если переход из начальной точки в конечную выполнен при числе переходов с одной траектории на другую больше одного, то в данной задаче управление выбрано не оптимальным.
Пусть конечным состоянием является начало координат (y10=0, y20=0). Полуветви парабол, по которым можно попасть в точку (y10=0, y20=0), образует линию переключения, которую обозначим MON. Уравнение линии переключения в этом случае имеет вид:
y1=-y22/2*sign(y2(t))(28)
Задача на максимальное быстродействие будет решена, если из начальной точки выбирается управление, при котором обеспечивается выход на линию переключения. При встрече с этой линией осуществляется переключение знака управления на противоположный, что обеспечивает переход на линию переключения и спуск по ней в конечную точку.
Причем, если начальная точка расположена правее линии переключения, то движение из нее нужно осуществлять при U=1. Если начальная точка расположена правее линии переключения, то движение из нее нужно осуществить при U=-1.
Таким образом, оптимальное управление формируется условиями:
·Если y1(t)<-y22/2*sign(y2(t)), то Uопт(t)=1(29)
·Если y1(t)>-y22/2*sign(y2(t)), то Uопт(t)=-1(30)
В том случае, если начальная точка расположена на линии переключения, т.е. y1(t)=-y22/2*sign(y2(t)), то Uопт(t)=1 (31); если y2(t)<0 или Uопт(t)=-1, если y2(t)>0 (32).
В этом случае переключения знака управления не требуется. Математически это означает, что за время перехода из начальной точки в конечную функция ?2(t) не изменяла знак.
Результатом моделирования оптимального перехода системы из заданного начального в заданное конечное состояние являются графики движения системы в фазовом пространстве и во временной области. Начальное и конечное состояния системы задаются преподавателем.
Просмотров: 31349
Найти в Wikkipedia статьи с фразой: Экономико-математическое моделирование