Оптимизация методом сканирования
Методы и программы одномерной оптимизации.
В инженерной и научной деятельности достаточно часто возникает необходимость в одномерной оптимизации, т.е. необходимость в нахождении экстремума заданной функции в некотором диапазоне значений аргумента.
Наиболее простым, в смысле понимания, является метод сканирования, состоящий в последовательном анализе значений исследуемой функции.
Далее приводится вариант программы на языке С (компилятор BORLAND C++ 3.1)
//****************************************************************
// Программа поиска ГЛОБАЛЬНОГО экстремума функции
// МЕТОДОМ СКАНИРОВАНИЯ С УТОЧНЕНИЕМ.
// Суть метода состоит в двукратном сканировании
// зоны интереса функции.
// - первый раз от начала до конца зоны интереса с
// заданным начальным шагом сканирования (D). Анализируется
// модуль значения функции. Так находится ГЛОБАЛЬНЫЙ
// ЭКСТРЕМУМ (максимум или минимум) - EXTR и XEXTR.
// - второй раз от XEXTR-D до XETR+D с шагом сканирования
// равным D/(заданный коэф. уменьшения шага сканирования).
// Программа при вводе -1 2 0.1 100 выдает
// 0.571000000000000 0.999999979262692
//****************************************************************
#include <conio.h>
#include <stdio.h>
#include <math.h>
float f(float x)
{
return sin(1+x);
}
int main(void)
{
double d,x,xk,extr=0,xextr,y;
int k,k1=1;
puts("Введите левую и правую границы зоны анализа, шаг движения и делитель шага ");
scanf("%lf%lf%lf%d",&x,&xk,&d,&k);
while (k1!=2*k)
{
(k1!=1)?(x=xextr-d,xk=xextr+d,d/=k1):(0);
x-=d;
while ((x+=d)<=xk)
((y=fabs(f(x)))>extr)?(extr=y, xextr=x):(0);
(k1==k)?(k1=2*k):(k1=k);
}
printf("%20.15lf %20.15lf",xextr,extr);
}
Далее приводится текст программы на языке PASCAL. Проверено на компиляторах
PASCALABC 3.1 и PASCALABC.NET 2.2.
//****************************************************************
// Программа поиска ГЛОБАЛЬНОГО экстремума функции
// МЕТОДОМ СКАНИРОВАНИЯ С УТОЧНЕНИЕМ.
// Суть метода состоит в двукратном сканировании
// зоны интереса функции.
// - первый раз от начала до конца зоны интереса с
// заданным начальным шагом сканирования (D). Анализируется
// модуль значения функции. Так находится ГЛОБАЛЬНЫЙ
// ЭКСТРЕМУМ (максимум или минимум) - EXTR и XEXTR.
// - второй раз от XEXTR-D до XETR+D с шагом сканирования
// равным D/(заданный коэф. уменьшения шага сканирования).
// Программа при вводе -1 2 0.1 100 программа выдает
// 0.571000000000000 0.999999979258613
//****************************************************************
var
d,x,xk,extr,xextr,y:real;
k,k1:integer;
function f(x:real):real;
begin
f:= sin(1+x);
end;
begin
extr:=0;
k1:=1;
writeln('Введите левую и правую границы зоны анализа, шаг движения и делитель шага');
readln(x,xk,d,k);
while (k1<>2*k) do
begin
if(k1<>1) then
begin
x:=xextr-d;
xk:=xextr+d;
d:=d/k1;
end;
while (x<=xk) do
begin
y:=abs(f(x));
if(y>extr) then
begin
extr:=y;
xextr:=x;
end;
x:=x+d;
end;
if(k1=k) then k1:=2*k else k1:=k;
end;
writeln(xextr:20:15,extr:20:15);
end.
Далее приводится текс программы на языке QBASIC в компиляторе
MS-DOS QBASIC 1.0.
DECLARE FUNCTION FUN# (x AS DOUBLE) '**************************************************************** ' Программа поиска ГЛОБАЛЬНОГО экстремума функции ' МЕТОДОМ СКАНИРОВАНИЯ С УТОЧНЕНИЕМ. ' Суть метода состоит в двукратном сканировании ' зоны интереса функции. ' - первый раз от начала до конца зоны интереса с ' заданным начальным шагом сканирования (D). Анализируется ' модуль значения функции. Так находится ГЛОБАЛЬНЫЙ ' ЭКСТРЕМУМ (максимум или минимум) - EXTR и XEXTR. ' - второй раз от XEXTR-D до XETR+D с шагом сканирования ' равным D/(заданный коэф. уменьшения шага сканирования). ' Программа при вводе -1,2,0.1,100 выдает ' 0.57100000000000 0.999999979258613 '**************************************************************** DIM d AS DOUBLE, x AS DOUBLE, xk AS DOUBLE, xextr AS DOUBLE DIM extr AS DOUBLE, y AS DOUBLE, k AS INTEGER, k1 AS INTEGER 'Тело программы CLS 'Очиcтка экрана extr = 0 k1 = 1 'Ввод исходных данных PRINT "Введите левую и правую границы зоны анализа, шаг движения и делитель шага "; PRINT INPUT x, xk, d, k WHILE (k1 <> 2 * k) IF (k1 <> 1) THEN x = xextr - d: xk = xextr + d: d = d / k1 WHILE (x <= xk) y = ABS(FUN(x)) IF (y > extr) THEN extr = y: xextr = x x = x + d WEND IF (k1 = k) THEN k1 = 2 * k ELSE k1 = k WEND PRINT USING "##.############## ##.###############"; xextr; extr END ' Функция собственно вычисления функции FUNCTION FUN# (x AS DOUBLE) FUN = SIN(1 + x) END FUNCTION
Оптимитизация методом половинного деления и золотого сечения
Ниже приведены на различных языках программы поиска выделенного экстремума по двум алгоритмам. Поясним суть использованных алгоритмов.
Половинное деление
Вычисляется значение анализируемой функции в двух точках, отстоящих от середины зоны локализации искомого экстремума вправо и влево на величину, равную половине заданной точности поиска положения экстремума. Затем сравниваются значения функции в этих точках. Из дальнейшего анализа исключается та половинка, значение функции в которой меньше. И так до достижения такого состояния, когда зона анализа станет меньше заданной точности (размера) зоны анализа. На каждом этапе уменьшения зоны анализа вдвое анализируемая функция вычисляется два раза.
Метод золотого сечения
В некотором смысле данный метод подобен предыдущему (подробное описание самой идеи золотого сечения, там он используется для поиска выделенного корня). Достоинством данной реализации метода поиска экстремума по сравнению с предыдущим является то, что на каждом этапе сокращения зоны анализа значение анализируемой функции вычисляется только один раз.
Ниже приведены тексты соответствующих программ на различных языках программирования.
Далее приводится вариант программы на языке С (компилятор BORLAND C++ 3.1)
//******************************************************************
// MAX02.CPP
// Программа поиска МАКСИМУМА функции методом ПОЛОВИННОГО ДЕЛЕНИЯ
// и ЗОЛОТОГО СЕЧЕНИЯ.
// Программа является переработанным вариантом программы MAX01.CPP.
// Достоинством этой реализации является то, что на каждом шаге
// значение функции вычисляется только ОДИН РАЗ в золотом сечении и
// два раза в половинном делении.
// Работа программы:
// Введите aa, bb: -1 2
// Точн._задан. ПОЛОВИННОЕ_ДЕЛЕНИЕ ЗОЛОТОЕ_СЕЧЕНИЕ
// XMAX вызов_функции XMAX вызов_функции
// 1.0e-04 0.5707923029 54 0.5708101855 24
// 1.0e-05 0.5707949120 67 0.5707965244 29
// 1.0e-06 0.5707962892 77 0.5707965244 33
// 1.0e-07 0.5707963235 88 0.5707963447 38
// 1.0e-08 0.5707963266 101 0.5707963285 43
// 1.0e-09 0.5707963269 112 0.5707963270 48
// 1.0e-10 0.5707963270 121 0.5707963276 51
//******************************************************************
#include <iostream.h>
#include <conio.h>
#include <math.h>
#include <stdio.h>
int k,k1;
double f(double x)
{
k++; k1++;
return sin(1 + x);
}
void main()
{
double a,aa, b, bb, c, d, e, f1, f2,dd;
clrscr();
cout<<"Введите aa, bb: ";
cin>>aa>>bb;
cout<<"Точн._задан. ПОЛОВИННОЕ_ДЕЛЕНИЕ ЗОЛОТОЕ_СЕЧЕНИЕn";
cout<<" XMAX вызов_функции XMAX вызов_функции ";
for(e=1e-4;e>1e-10;e/=10)
{
//**********************
// ПОЛОВИННОЕ ДЕЛЕНИЕ
//**********************
k=0;
a=aa;
b=bb;
while(fabs(a-b)>e)
((f1=f((d=(a+b)/2)-e/2))>(f2=f(d+e/2)))? b=d:a=d;
//12 34 5 5 4 32 6 7 761
dd=d;
//****************************
// ЗОЛОТОЕ СЕЧЕНИЕ
//*****************************
k1=0;
a=aa;
b=bb;
f1=f(c=a+(b-a)*(sqrt(5.)-1)/2);
f2=f(d=a+b-c);
while((fabs(a-b)>e) && (fabs(f1-f2)>e))
(f1<f2)?(f1=f2,f2=f(d=a+(b=c)-(c=d))):(f2=f1,f1=f(c=b+(a=d)-(d=c)));
//1 1 2 3 4 4 5 532 6 7 8 8 9 976
printf("n%5.1e %15.10f %6d %25.10f %5d",e,dd,k,c,k1);
}
getch();
}
Далее приводится текс программы на языке QBASIC в компиляторе
MS-DOS QBASIC 1.0.
DECLARE FUNCTION FUN# (x AS DOUBLE) '****************************************************************** ' MAX.CPP '****************************************************************** ' MAX02.CPP ' Программа поиска МАКСИМУМА функции методом ПОЛОВИННОГО ДЕЛЕНИЯ ' и ЗОЛОТОГО СЕЧЕНИЯ. При разных значениях точности. ' Достоинством этой реализации является то, что на каждом шаге ' значение функции вычисляется только ОДИН РАЗ в золотом сечении и ' два раза в половинном делении. ' Работа программы: ' Введите границы зоны aa, bb: -1,2 'Точн._задан. ПОЛОВИННОЕ_ДЕЛЕНИЕ ЗОЛОТОЕ_СЕЧЕНИЕ ' XMAX вызов_функции XMAX вызов_функции '0.1D-03 0.5707702637 39 0.5835921350 9 '0.1D-04 0.5707988739 50 0.5742752750 12 '0.1D-05 0.5707967281 59 0.5720758626 15 '0.1D-06 0.5707962811 67 0.5712357619 17 '0.1D-07 0.5707963314 77 0.5709148720 19 '0.1D-08 0.5707963614 86 0.5708391201 22 '0.1D-09 0.5707967151 94 0.5708101855 24 '0.1D-10 0.5708049759 104 0.5707991335 26 '****************************************************************** '********************************************************************* DIM a AS DOUBLE, aa AS DOUBLE, b AS DOUBLE, bb AS DOUBLE DIM c AS DOUBLE, d AS DOUBLE, e AS DOUBLE, f1 AS DOUBLE, i AS INTEGER DIM f2 AS DOUBLE, dd AS DOUBLE DIM SHARED k AS INTEGER, k1 AS INTEGER CLS INPUT "Введите границы зоны анализа aa и bb ", aa, bb PRINT "Точн._задан. ПОЛОВИННОЕ_ДЕЛЕНИЕ ЗОЛОТОЕ_СЕЧЕНИЕ" PRINT " XMAX вызов_функции XMAX вызов_функции " PRINT e = .001 ' aa = -1 ' bb = 2 FOR i = 1 TO 8 e = e / 10 '********************** ' Половинное деление '********************** f1 = 1 f2 = 0 k = 0 a = aa b = bb DO WHILE (ABS(a - b) > e) d = (a + b) / 2 IF (FUN(d - e / 2) > FUN(d + e / 2)) THEN b = d ELSE a = d LOOP dd = d '**************************** ' Золотое сечение '***************************** k1 = 0 a = aa b = bb c = a + (b - a) * (SQR(5#) - 1) / 2 f1 = FUN(c) d = a + b - c f2 = FUN(d) WHILE ((ABS(a - b) > e) AND (ABS(f1 - f2) > e)) IF (f1 < f2) THEN f1 = f2: b = c: c = d: d = a + b - c: f2 = FUN(d) ELSE f2 = f1: a = d: d = c: c = b + a - d: f1 = FUN(c) END IF WEND PRINT USING "#.#^^^^ ##.########## #### ##.########## ####"; e; dd; k; c; k1 NEXT i END ' Собственно вычисление функции FUNCTION FUN# (x AS DOUBLE) k = k + 1 k1 = k1 + 1 FUN = SIN(1 + x) END FUNCTION
Ниже приводится текст этой же программы составленный для интерпретатора Python 3.4.3.
#******************************************************************
# MAX1.PY
#******************************************************************
# Программа поиска МАКСИМУМА функции методом ПОЛОВИННОГО ДЕЛЕНИЯ
# и ЗОЛОТОГО СЕЧЕНИЯ. При разных значениях точности.
# Достоинством этой реализации является то, что на каждом шаге
# значение функции вычисляется только ОДИН РАЗ в золотом сечении и
# два раза в половинном делении.
# Работа программы:
# Введите границы зоны aa, bb: -1,2
#Точн._задан. ПОЛОВИННОЕ_ДЕЛЕНИЕ ЗОЛОТОЕ_СЕЧЕНИЕ
# XMAX вызов_функции XMAX вызов_функции
#1.0e-04 0.5707702637 39 0.5835921350 9
#1.0e-05 0.5707988739 50 0.5742752750 12
#1.0e-06 0.5707967281 59 0.5720758626 15
#1.0e-07 0.5707962811 67 0.5712357619 17
#1.0e-08 0.5707963314 77 0.5709148720 19
#1.0e-09 0.5707963614 86 0.5708391201 22
#1.0e-10 0.5707967151 94 0.5708101855 24
#1.0e-11 0.5708049759 104 0.5707991335 26
#******************************************************************
#*********************************************************************
#import pdb; pdb.set_trace()
import math
# Собственно вычисление функции
def F(x):
global k,k1
k+=1
k1+=1
return math.sin(x+1)
(aa,bb)= list(map(float, input('Введите границы зоны анализа aa и bb ').split()))
#aa=-1
#bb=2
print('Точн._задан. ПОЛОВИННОЕ_ДЕЛЕНИЕ ЗОЛОТОЕ_СЕЧЕНИЕ')
print(' XMAX вызов_функции XMAX вызов_функции ')
e=.001
k=0
k1=0
for i in range(1,9):
e/=10
#**********************
# Половинное деление
#**********************
k=0
a=aa
b=bb
while (abs(a-b)>e):
d=(a+b)/2
if (F(d-e/2)>F(d+e/2)): b=d
else: a=d
dd=d
#****************************
# Золотое сечение
#*****************************
k1=0
a=aa
b=bb
c=a+(b-a)*(math.sqrt(5.)-1)/2
f1=F(c)
d=a+b-c
f2=F(d)
while ((abs(a-b)>e) and (abs(f1-f2)>e)):
if(f1<f2): f1=f2; b=c; c=d; d=a+b-c; f2=F(d)
else: f2=f1; a=d; d=c; c=b+a-d; f1=F(c)
text='%5.1e %12.10f %4d %12.10f %4d'%(e,dd,k,c,k1)
print(text)