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

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ

Государственное образовательное учреждение

высшего профессионального образования

«МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

ПИЩЕВЫХ ПРОИЗВОДСТВ»

         Кафедра «Автоматика и Электротехника».

Курсовая работа.

«Численные методы решения задач управления технологическими процессами». 



Группа: 07-ИУ-4

Студент: Коняхин Е. И.

Преподаватель: Михайлов А. В.











Москва 2010.



Поиск глобального максимума методом равномерного поиска.

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



Пусть  заданная погрешность определения максимального значения  .

Тогда для реализации алгоритма поиска следует определить значение  в



точках, равномерно относящихся друг от друга на расстоянии  , т.е. в точках:



Из полученных значений показателя качества  выбирается наибольшее значение (глобальный максимум).







началоАлгоритм расчета.



f=

x

a,b,h













 





x=x+h





	

конец



f=f(x)

	ДА









		

	

	

	

	Результаты расчета.

	Целевая функция имеет вид :  

	Интервал (-100;100)

	Шаг H

	Значение X

	Значение F

	Кол-во вычислений

	1

	1

	23

	21

	0,1

	0,19999998

	23

	201

	0,01

	0,959999

	23,056

	2001

	

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

	

	Одномерная оптимизация методом дихотомии.

	Этот метод используется для поиска экстремума класса унимодальных функций. Идея метода проста – делить интервал [a,b] , где расположена точка экстремума  , пополам

	

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

	  

	 

	Таким образом, в результате двух вычислений  , промежуток , где содержится экстремум сократится в двое. Следующая пара измерений производится в районе серидины нового интервала неопределенности. Вычисления производятся до тех пор, пока на k-м шаге, после 2k вычислений  длина интервала неопределенности , где находится оптимум, не станет меньше или равна .

	

	

	

	

	

	

	

	

	

	

	

	Алгоритм расчёта.

	



x1=

x2=

f1=f(x1)

f2=f(x2)

а=x1

b=x2

f2<f1

конец







начало

a,b,h

d=

	















		ДА









		

	       ДА





		

	Результаты расчета.

	Целевая функция имеет вид :  

	Интервал (-100;100)

	Погрешность Е

	Значение Х

	Значение F

	Кол-во итераций

	Кол-во вычислен.

	0.1

	1,06640625

	22,600226978

	14

	7

	0.1

	0,954638671

	23,054333993

	20

	10

	0.01

	0,959232788

	23,05589651

	26

	13

	0.001

	0,96173743

	23,056110564

	34

	17

	

	             Вывод : как видно метод  дихотомии позволяет довольно быстро попадать в район оптимума. И требует меньшего числа расчётов по сравнению с некоторыми другими методами (например : методом равномерного приближения).

	

	

	

	

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

	          Интервал неопределенности делится на три отрезка, причем внутренние точки располагаются симметрично по отношению к крайним .

	          Берутся пробные точки и располагаются следующим образом : 

	

	

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

	Поиск оптимума завершается, если после k- го шага длина интервала неопределенности станет меньше или равна .

	

	Алгоритм расчёта.

	

начало

	

x1=b-0.62(b-a)

f1=



b=x2

x2=x1

f2=f1

конец







(b-a)≤ɛ



=



a=x2

x1=x2

f1=f2



f1<f2

x1=a+0.62(b-a)

f1=



x1=b-0.62(b-a)

f1=

a,b,h















		нет

		нет





















Результаты расчета.

	Целевая функция имеет вид :  

Интервал (-100;100)

Погрешность Е

Значение Х

Значение F

Кол-во итераций

Кол-во вычислен.

1

0,895427987

22,910241869

9

8

0.1

0,94547427929

23,046790768

14

13

0.01

0,96349811902

23,055988462

18

17

0.001

0,96148191718

23,05610953

23

22



         Вывод : преимуществом этого метода над методом дихотомии является то, что на каждом шаге вычисляется лишь одно значение , а не два.



Одномерная оптимизация методом поразрядного приближения.

         Метод обладает высоким быстродействием. Это достигается тем, что используется алгоритм с переменным шагом поиска. Задаем интервал [a,b] , содержащий внутри себя точки максимума :  

         Задается начальное значение  и вычисляется . Задается начальный шаг поиска h и кратность изменения шага k в районе оптимума. Производится поиск максимума. Поиск из начальной точки x= осуществляется с постоянным шагом h , после каждого шага вычисляется значение критерия  , оно сравнивается с предыдущим и в случае улучшения критерия шаги продолжаются. Движение к оптимуму с неизменным шагом h продолжается до тех пор, пока очередной шаг не  окажется неудачным. После этого поиск максима продолжается из последней точки в обратном направлении с шагом в k раз меньше прежнего. Эта процедура будет продолжаться до тех пор, пока не выполнится условие: 

          , где - заданная погрешность определения оптимума.































началоАлгоритм расчета.

  



=f

x=x+h

f=f(x)

h=-h/k

конец









 ,h, , k



















да



		нет













		



Результаты расчета.

	Целевая функция имеет вид :  

	X=-100;h=10;k=10

	Погрешность Е

	Оптимальная точка х

	Оптимальное значение ф-ции

	Кол-во итераций

	Кол-во вычислений

	1

	0,895427987

	22,910241869

	2

	23

	0.1

	0,94547427929

	23,046790768

	3

	35

	0.01

	0,96349811902

	23,055988462

	4

	51

	0.001

	0,96148191718

	23,05610953

	5

	65

	

	

	

	X=-100;h=10;e=0.01

	Кратность k

	Оптимальная точка х

	Оптимальное значение ф-ции

	Кол-во итераций

	Кол-во вычислений

	2

	0,95703125

	23,0553

	11

	51

	5

	0,9632

	23,0560

	6

	46

	10

	0,96

	23,0560

	4

	51

	20

	0,96125

	23,0560

	4

	65

	50

	0,96

	23,0560

	3

	101

	

	X=-100;k=10;e=0.01

	Шаг h

	Оптимальная точка х

	Оптимальное значение ф-ции

	Кол-во итераций

	Кол-во вычислений

	0.1

	0,95999

	23,05601

	2

	1028

	0.5

	0,96

	23,05601

	3

	231

	1

	0,96

	23,05601

	3

	123

	5

	0,96

	23,05601

	4

	53

	10

	0,96

	23,05601

	4

	51

	50

	0,96

	23,05601

	5

	57

	100

	0,96

	23,05601

	5

	48

	

	          Вывод: метод является эффективным для измерения оптимума унимодальной функции, причем изменение шага поиска или кратности уменьшения шага ( при неизменной погрешности вычисления на результат практически не влияет).

Одномерная оптимизация методом квадратичной интерполяции.

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

(х) = ах2 + Ьх + с

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

Экстремум функции  fап (х) как известно расположен в точке: = -Ь/2а. 

Положим, что окрестность некоторой исходной точки х=х1 области определения f0(х) аппроксимирована полиномом fап (х). Задача поиска заключается в определении смещения

 = х°ап – х1

Которое приводит из исходного состояния х = х1, ближе к экстремуму х = х°. Если f0(х)  строго квадратичная функция, то смещение  после первого шага сразу приведет к. В противном случае достижение х° требует выполнения итерационной процедуры. Для определения смещения  нужно определить коэффициенты параболы. Для этого необходимо вычислить значение f0(х) в трех точках. Пусть вычисление производится в исходном состоянии х = х1 и в точках, ,  и при этом получено три значения этой функции

,

 где h - полуинтервал интерполяции, малая постоянная величина. Подставляя эти значения в уравнение (х), получаем систему из трех линейных уравнений с тремя неизвестными а, Ь, с:

а(х1 - h)2 + Ь(х1 - h) + с =а(х1 - h)2 + Ь(х1 - h) + с =

а*х12 + Ь*х1 + с =

а(х1 + h)2 + Ь(х1 + h) + с =

Для того, чтобы система имела решение, необходимо чтобы ее определитель не был равен нулю. Это условие выполняется, так как определитель равен: = - 2h30 так как . Решая систему уравнений, получаем интересующие нас значения параметров а, Ь, с подставляя их в формулу  находим положение экстремума параболы

х°ап= х1 + h(-)/2(- 2+)

Зная коэффициенты а, Ь, с можно определить и экстремальное значение функции по формуле, которая является оценкой экстремума критерия (х).



Теперь следует проверить, действительно ли найден экстремум. Для этого достаточно вычислить значение функции цели (х) в предполагаемом экстремуме х=х1+Δх - х°ап и сопоставить его с оценкой. Если эти величины отличаются не более чем на ɛ т. е:

|( х°ап )-(х°ап )|

, где ɛ заданная погрешность определения экстремума. При этом  = х1. Если условие не выполняется, тогда следует процесс поиска; т.е. выполнить следующий цикл, но уже построение

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

Алгоритм расчета.











Результаты расчета.

	Целевая функция имеет вид :  

Нач. знач. X=-100,H=0.5

Погрешность Е

Значение Х

Значение F

Кол-во итераций

Кол-во вычислений

1

(-)2,19360741

(-)919,076558

10

30

0.1

0,8912446

22,8921666

14

45

0.01

0,79728604

22,27161267

16

48

0.001

0,7960595

22,2612358

17

51

Нач. знач. Х=-100, Е=0.1

Шаг Н

Знач Х

Знач F

Кол-во итераций

Кол-во вычислений

Увеличение шага

0,3

0,901465

22,93463

25

75

0,5

0,797286

22,27161

16

48

0,8

0,6115913

20,33949

35

105

Уменьшение шага

0.5

0,79728604

22,27161267

16

48

0.4

0,8540667

22,69232

20

60

0.3

0,901465

22,934634

25

75

0.2

0,936694

23,034198

41

123

0.1

0,961479

23,056109

31

93

0.02

0,961661

23,05611

25

75

                

          Вывод: расчеты показали, что изменение погрешности определения экстремума ɛ, практически не влияет на точность вычисления в то время, как изменение шага поиска h оказывает значительное влияние. При уменьшении шага точность вычислений улучшается и наоборот, при увеличении шага уменьшается. И в конечном итоге, когда шаг поиска слишком велик для того, чтобы с помощью итерационной процедуры уточнения значений получить результат с заданной погрешностью, программа отказывается производить вычисления.

Оптимизация методом наискорейшего спуска.

Метод наискорейшего спуска предназначен для поиска минимума. Данный метод отличается от метода градиента правилом определения коэффициента шага. Сначала выделяется начальная точка. В пространстве X могут быть выделены области притяжения каждого из локальных минимумов.

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













Алгоритм расчета.



начало





n,,h,ɛ,l



	

	

	

	

	

	

	

	Результаты расчета.

	Целевая функция имеет вид :  

Н=1,E=0.01

Приращение L

Оптим. зн. XI

Оптим. зн. Х2

Оптим. зн. ХЗ

Оптим. зн. F

0,5

-10,75

17,25

-12,75

0,1875

0.1

-10,95

17,04

-12,95

0,0074

0.01

-10,995

17,005

-12,994

7,5000001E-7

0.001

-10,99

17

-12,99

7,49752E-9

	

L=0.0001, Е=0.1

Шаг h

Оптим. зн. XI

Оптим. зн. Х2

Оптим. зн. ХЗ

Оптим. зн. F

10

-100

100

-100

31979

5

-100

100

-100

31979

1

-10,9

17

-12,9

7,5

0,5

-10,9

17

-12,9

7,49934E-9

0,1

-10,9

17

-12,9

7,49934E-9

	

Н=0.5, L=0.0001

Погрешность Е

Оптим. зн. X 1

Оптим. зн. Х2

Оптим. зн. ХЗ

Оптим. зн. F

1

-121,99

65,9

-125,9

31978,93

0.1

-10,99

17

-12,99

7,49752Е-9

0.01

-10,99

17

-12,99

7,49752Е-9

0.001

-10,99

17

-12,99

7.49752Е-9

	

Оптимизация методом линейного программирования.

	f0(x)=4x+3y

	

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

1)3x-y≤4

2)x-2y≤-7

3)3x+y≤21

4)-x+4y≤6

Следовательно искомые ограничения:

1)-3x+y-4

2)-x+2y≤7

3)3x+y≤21

4)x-4y≤-6

Точки min и max:

Amin (2;2)   Cmax (5;6)

min и max функции:

	f0(x)min=14

	f0(x)max=38

Расчет производится в приложении МАТLАВ.

»f=[4,3]

 f = 4     3

»A=[-3,l;-l,2;3,l;l,-4] 

A =

-3     1

-1     2

3     1

1    -4

» B=[-4;7;21;-6]

 B =

-4

7

21

-6

» [x,y,z]=linprog(f,A,B) Optimization terminated, 

x =

2.0000

2.0000

y =

14.0000

z = 1

Решение задачи нелинейного программирования.

	Задача №2.

	Имеются три продукта n1,n2,n3 разной цены. Каждый из них содержит определенное количество питательных ингредиентов, причем для нормального потребления требуется u1≥250; u2≥60; u3≥100; u4≥220. 

	Минимизировать затраты на приобретение продукта.

	

	N1

	N2

	N3

	U1

	4

	6

	15

	U2

	2

	2

	0

	U3

	5

	3

	4

	U4

	7

	3

	12

	Цена за единицу

	44

	35

	100

	

Расчет производится в приложении МАТLАВ.

a=[-4 -6 -15;-2 -2 0;-5 -3 -4;-7 -3 -12]; »b=[-250-60-100-220]; 

»f=[44 35 100] 

f =

44    35   100 

» [x,y,z]=linprog(f,a,b) Optimization terminated. x= 13.2143

16.7857

6.4286

y =

1.8118e+003 z= 1

	

	Задача №3 Вариант№13



 N1

 N2

  N3

  N4

 N5

В1

А1

2

5

7

9

12

4000

А2

4

9

16

21

28

6000

АЗ

7

12

18

27

34

8000

А4

10

18

26

34

42

10000

А5

14

23

32

44

51

12000

А6

18

29

39

50

63

14000

С

10

11

12

13

14



	

Расчет производится в приложении МАТLАВ.

	

a=[2,5,7,9,12;4,9,16,21,28;7,12/18,27,34;10,18,26,34,42;14l23,32,44,51;18,29,39,50,63]

b=[4000;60000;8000;10000;12000;14000]

lb=zeros(5,l)

f=[-10-ll -12-13-14]

»[x,y,z]=linprog(f,a,b,[],[],lb)

Optimization terminated.

x =

777.7778 

0.0000 

0.0000 

0.0000 

0.0000

y =

-7.7778e+003

z =

1



	

Задача №4 Вариант №8

12.

Уваривание рецептурной смеси в змеевиковом уварочном аппарате.

Общий вид зависимости:  τвар=f(P пара)

Результаты моделирования: y=178-865,5*x+1629,4*x^2-1025*x^3

Минимизировать время уваривания смеси до влажности 16% при ограничении на давлении в магистрали пара, обусловленного экономической эффективность.

Ограничения: 0,2 <Р пара<0,7 Мпа; f`(Р пара)≤-44

Расчет производится в приложении МАТLАВ.

function f=myfun(x)

function[c,ceq]=ogr(x); c=-865,5+3258,8*x-3075*x^2+44; 

ceq=[];

»[x,y,z]=fmincon(@myfun,[1],[],[],[],[],[0,7],[0,2],@ogr)

x =

0

Y =

-687

z =

1

	

	13.

Уваривание рецептурной смеси в змеевиковом уварочном аппарате.

Общий вид зависимости:  y отк=f(P апп)

Результаты моделирования: y=47,5-2,01*x+0,03*x^2-0,000015*x^3

Повысить качество смеси, выраженное в заданной взякости 20 Па*с в зависимости от ограничения для температуры в аппарате и частотой вращения мешалки. 

Ограничения: 60 <Т апп<76 °С; 15<w лоп<35 об/мин;



Расчет производится в приложении МАТLАВ.

function f=myfun(x)

f=47,5-2,01*x+0,03*x^2-0,00015*x^3;

function[c,ceq]=ogr(x); c=-2,01+0,06*x-0,00045*x^2; ceq=[];

»[x,y;z]=fmincon(@myfun,[l],[],[],[],[],[65],[93],@ogr)

x =

65.0001

y =

47

z =

1