Экономико-математическое моделирование

ЛАБОРАТОРНАЯ РАБОТА №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 Пример выполнения задания


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

Исследовать на экстремум функции:

  1. y=10-48x+16x2;
  2. 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)


  1. Необходимое условие экстремума:
  1. Допустимое условие экстремума:

Так как порядок, впервые не обращающейся в нуль производной равен 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) не изменяла знак.

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


Теги: Экономико-математическое моделирование  Практическое задание  Менеджмент
Просмотров: 31256
Найти в Wikkipedia статьи с фразой: Экономико-математическое моделирование
Назад