close

Вход

Забыли?

вход по аккаунту

?

Курсовик (отчет)

код для вставкиСкачать
 Министерство образования и науки РФ
Государственное образовательное учреждение высшего профессионального образования "Ивановский государственный энергетический университет им. В. И. Ленина"
Кафедра высокопроизводительных систем
Отчет по курсовой работе
по курсу "Методы оптимизации"
"Методы безусловной оптимизации второго порядка: методы Ньютона, Ньютона-Рафсона и Марквардта"
Выполнил: студент гр. 3-41х
Громов А.А
Проверил: Филатов Е.Ю.
Иваново 2011
Содержание
Задание.........................................................................................................3
Теоретическая часть......................................................................................3
1. Описание методов............................................................................3
2. Описание алгоритмов на псевдокоде.....................................................6
Практическая часть ......................................................................................8
1. Описание интерфейса программы..... ....................................................8
2. Используемые в программе классы, переменные, функции и т.п...................8
3. Результаты работы программы (рисунки, таблицы, ответы) ........................ 10
Библиографический список ............................................................................................ 13
Приложение 1. Код программы............................................................................................. 14
Задание
Написать программу, которая реализует методы Ньютона, Ньютона-Рафсона и Марквардта для поиска минимума функции Начальное приближение, точность решения и метод оптимизации выбираются пользователем. Программа должна строить график линий уровня функции f(x,y), указывать на нем промежуточные точки и найденную точку минимума. Если точка минимума не была найдена, необходимо выдать соответствующее сообщение.
Теоретическая часть
Описание методов
Метод Ньютона
Стратегия метода Ньютона состоит в построении последовательности точек ,
k=0, 1,... таких, что , k=0, 1,...
Точки последовательности вычисляются по правилу:
где - задается пользователем, а направление спуска вычисляется по формуле:
Такой выбор направления спуска гарантирует выполнение требования , при условии, что матрица Гессе положительно определена. Чтобы обеспечить выполнение требования даже в тех случаях, когда для каких либо значений матрица Гессе не окажется положительно определенной, рекомендуется для соответствующих значений k вычислить точку по методу градиентного спуска с выбором величины шага из условия К методу Ньютона можно придти из следующих соображений.
Необходимым условием экстремума функции многих переменных f(x) в точке является равенство нулю ее градиента в этой точке:
Функция f(х) в окрестности точки может быть разложена в ряд Тейлора с точностью до членов второго порядка:
Направление определяется из необходимого условия экстремума первого порядка:
.
Имеем При выполнении условия о положительной определенности матрицы получим:
Таким образом, решение задачи методом Ньютона предполагает построение последовательности минимумов аппроксимирующих квадратичных функций .
В случае двух переменных данный метод имеет следующую геометрическую интерпретацию.
Построение последовательности заканчивается в точке, для которой Если функция f(x) является квадратичной, то, независимо от начального приближения и степени овражности, с помощью метода Ньютона ее минимум находится за один шаг. Действительно, в этом случае аппроксимирующая функция совпадает с f(x) и метод Ньютона при положительной определенности матрицы Гессе отыщет минимум за один шаг.
Метод Ньютона-Рафсона
Существенным недостатком метода Ньютона является зависимость сходимости для невыпуклых функций от начального приближения . Если находится достаточно далеко от точки минимума, то метод может расходиться, т. е. при проведении итерации каждая следующая точка будет более удаленной от точки минимума, чем предыдущая. Для ослабления влияния начальной точки на сходимость метода используется метод Ньютона-Рафсона. Стратегия метода Ньютона-Рафсона состоит в построении последовательности точек , k=0, 1, ... таких, что , k=0, 1, ...
Точки последовательности вычисляются по правилу:
где - задается пользователем, а направление спуска вычисляется по формуле:
Величина шага выбирается из условия , как в методе наискорейшего спуска.
Метод Марквардта
При использовании метода Ньютона или метода Ньютона-Рафсона вследствие накопления ошибок в процессе счета матрица Гессе на некоторой итерации может оказаться отрицательно определенной или ее нельзя будет обратить.
Для обеспечения положительной определенности матрицы Гессе в точках последовательности используется метод Марквардта (метод Ньютона с регулировкой матрицы). В методе Марквардта точки строятся по закону: где - последовательность положительных чисел, обеспечивающих положительную определенность матрицы . Обычно в качестве берется значение на порядок больше наибольшего элемента матрицы . Так в ряде стандартных программ полагается . Если
то , в противном случае Очевидно, что если , то метод Марквардта представляет собой метод Ньютона, а если велико, то поскольку при больших метод Марквардта близок к градиентному методу. Поэтому, подбирая значения параметра , можно добиться, чтобы метод Марквардта сходился.
Описание алгоритмов на псевдокоде
Метод Ньютона
Начальный этап
Задать ,.
Найти градиент функции в произвольной точке и матрицу Гессе H(x)
Положить k=0.
Основной этап
Шаг 1. Вычислить Шаг 2. Проверить выполнение критерия останова а) если критерий выполнен, расчет окончен, б) если критерий не выполнен, то перейти к шагу 3.
Шаг 3. Вычислить матрицу Гессе Шаг 4. Найти обратную матрицу Шаг 5. Проверить положительную определенность матрицы Шаг 6. Найти Шаг 7. Положить k= k +1 и перейти к шагу 1.
Метод Ньютона-Рафсона
Начальный этап
Задать ,.
Найти градиент функции в произвольной точке и матрицу Гессе H(x)
Положить k=0.
Основной этап
Шаг 1. Вычислить Шаг 2. Проверить выполнение критерия останова а) если критерий выполнен, расчет окончен, б) если критерий не выполнен, то перейти к шагу 3.
Шаг 3. Вычислить матрицу Гессе Шаг 4. Найти обратную матрицу Шаг 5. Проверить положительную определенность матрицы Шаг 6. Найти Шаг 7. Найти величину шага из условия Шаг 8. Вычислить Шаг 9. Положить k= k +1 и перейти к шагу 1
Метод Марквардта
Начальный этап
Задать ,.
Найти градиент функции в произвольной точке и матрицу Гессе H(x)
Положить .
Основной этап
Шаг 1. Вычислить Шаг 2. Проверить выполнение критерия останова а) если критерий выполнен, расчет окончен, б) если критерий не выполнен, то перейти к шагу 3.
Шаг 3. Вычислить матрицу Гессе Шаг 4. Вычислить Шаг 5. Найти обратную матрицу Шаг 6. Вычислить Шаг 7. Найти Шаг 8. Проверить выполнение условия :
а) если неравенство выполнено, то положить , k= k +1 и перейти к шагу 1;
б) если нет, то положить и перейти к шагу 4.
Практическая часть
Описание интерфейса программы
Интерфейс программы представляет собой приложение WindowsForms. Пользователь при помощи клавиатуры может изменить начальное приближение и точность решения. Выбор метода осуществляется при помощи кнопок. При нажатии на одну из них происходят вычисления соответствующим методом. Результат выводится в поля: x= и y= Так же строится график линий уровня с указанием начальной, промежуточных и конечной точек. В программе есть возможность масштабирования графика при помощи мыши. Если точка минимума не найдена, программа выдает соответствующее сообщение.
Используемые в программе классы, переменные, функции и т.п.
Основные функции
public double f(double x, double y) //функция вычисляющая значение целевой функции
public double dfdx(double x, double y) //функция вычисляющая значение первой производной целевой функции по х public double dfdy(double x, double y) //функция вычисляющая значение первой производной целевой функции по у
public double dfdx2(double x, double y) //функция вычисляющая значение второй производной целевой функции по х public double dfdy2(double x, double y) //функция вычисляющая значение второй производной целевой функции по у
public double dfdydx(double x, double y) //функция вычисляющая значение частной производной целевой функции сначала по у, потом по х
public double dfdxdy(double x, double y) //функция вычисляющая значение частной производной целевой функции сначала по х, потом по у
public double func(double t) //функция вычисляющая значение функции для определения шага в методе Ньютона-Рафсона
public double fib(double x1, double x2, int N) //функция вычисляющая минимум функции методом Фибоначчи
public void Newton(double x0, double y0,double ee) //функция вычисляющая минимум функции двух переменных методом Ньютона
public void Newton_Raphson(double x0, double y0,double ee,double a, double b) //функция вычисляющая минимум функции двух переменных методом Ньютона-Рафсона
public void Marquardt(double x0, double y0, double ee) //функция вычисляющая минимум функции двух переменных методом Марквардта
public void Graph() //функция для построения графика
Основные переменные
public PointPairList list1; //список для построения линий уровня
public PointPairList list2; //список для построения начальной, промежуточных и конечной точек
double[] grad; //массив для хранения текущего значения градиента
double[] d; //массив для хранения текущего направления спуска
double[] x и double[] y; //массивы для хранения предыдущего и текущего значения приближения по соответствующим координатам
double [,] Hisse; //двумерный массив для хранения матрицы Гессе
double Det; //переменная для хранения определителя матрицы
double m; //число, используемое в методе Марквардта
double t; //переменная для хранения значения шага
Результаты работы программы
Метод Ньютона
Метод Ньютона-Рафсона
Метод Марквардта
Пример неудачного вычисления точки минимума
Библиографический список
1. Математическое программирование в примерах и задачах/ Акулич И.Л. - М.: Высш. шк.,1986. 2. Методы оптимизации в примерах и задачах: Учеб. пособие / А.В.Пантелеев, Т.А.Летова. - М.: Высш. шк.,2002. 3. Методы оптимизации //Интернет-ресурс: http://www.theweman.info/topics/t4r10part1.html
Приложение 1. Код программы
public double f(double x, double y)
{
return x*Math.Exp(-Math.Pow(x,2)-Math.Pow(y,2));
}
public double dfdx(double x, double y)
{
return Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2)) - 2 * Math.Pow(x, 2) * Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2));
}
public double dfdy(double x, double y)
{
return -2 * x * y * Math.Exp (-Math.Pow(x, 2) - Math.Pow(y, 2));
}
public double dfdx2(double x, double y)
{
return -2 * x * Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2)) - 2 * (2 * x * Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2)) - 2*Math.Pow(x,3)*Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2)));
}
public double dfdy2(double x, double y)
{
return -2 *x *(Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2))-2 * Math.Pow(y, 2) * Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2)));
}
public double dfdydx(double x, double y)
{
return -2 * y * Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2)) + 4 * x * x * y * Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2));
}
public double dfdxdy(double x, double y)
{
return -2 * y * Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2)) + 4 * x * x * y * Math.Exp(-Math.Pow(x, 2) - Math.Pow(y, 2));
}
public double func(double t)
{
return (glx + t * gld1) * Math.Exp(-Math.Pow((glx + t * gld1), 2) - Math.Pow((gly + t * gld2), 2));
}
public double fib(double x1, double x2, int N)
{
int fa = 1, fb = 1, fbuf;
int Nbuf = 1;
while (Nbuf < N)
{
fbuf = fa + fb;
fa = fb;
fb = fbuf;
Nbuf++;
}
while (fb != 1)
{
double bb = (x1 + Convert.ToDouble(fa) / Convert.ToDouble(fb) * (x2 - x1));
double aa = x1 + x2 - bb;
if (func(aa) < func(bb))
x2 = bb;
else x1 = aa;
fa = fb - fa;
fb = fb - fa;
}
return (x1 + x2) / 2;
}
public void Newton_Raphson(double x0, double y0,double ee,double a, double b)
{
double t;
double[] grad = new double[2];
double[] d = new double[2];
double[] x = new double[2];
double[] y = new double[2];
double [,] Hisse=new double[2,2];
double Det, temp;
int flag;
bool CurFlag=false, PrevFlag=false, TrueResult=false;
x[0] = x0;
y[0] = y0;
list2.Add(x[0], y[0]);
int k = 0;
int counter = 0;
grad[0] = dfdx(x[k], y[k]);
grad[1] = dfdy(x[k], y[k]);
while (Math.Sqrt(Math.Pow(grad[0], 2) + Math.Pow(grad[1], 2)) > ee)
{
Hisse[0, 0] = dfdx2(x[k], y[k]);
Hisse[0, 1] = dfdxdy(x[k], y[k]);
Hisse[1, 0] = dfdydx(x[k], y[k]);
Hisse[1, 1] = dfdy2(x[k], y[k]);
Det = Hisse[0, 0] * Hisse[1, 1] - Hisse[0, 1] * Hisse[1, 0];
temp = Hisse[0, 0];
Hisse[0, 0] = Hisse[1, 1]/Det;
Hisse[1, 1] = temp/Det;
Hisse[0, 1] = -Hisse[0, 1]/Det;
Hisse[1, 0] = -Hisse[1, 0]/Det;
if (Hisse[0, 0] > 0 && Det > 0)
{
d[0] = -(Hisse[0, 0] * grad[0] + Hisse[0, 1] * grad[1]);
d[1] = -(Hisse[1, 0] * grad[0] + Hisse[1, 1] * grad[1]);
flag = 0;
}
else
{
d[0] = -grad[0];
d[1] = -grad[1];
flag = 1;
}
glx = x[k];
gly = y[k];
gld1 = d[0];
gld2 = d[1];
t = fib(a, b, 25);
x[k + 1] = x[k] + t * d[0];
y[k + 1] = y[k] + t * d[1];
list2.Add(x[k+1 ], y[k+1 ]);
if (counter == 0)
{
if ((Math.Sqrt(Math.Pow(x[k + 1] - x[k], 2) + Math.Pow(y[k + 1] - y[k], 2)) < ee) && (Math.Abs(f(x[k + 1], y[k + 1]) - Math.Abs(f(x[k], y[k]))) < ee))
PrevFlag = true;
}
else
{
if ((Math.Sqrt(Math.Pow(x[k + 1] - x[k], 2) + Math.Pow(y[k + 1] - y[k], 2)) < ee) && (Math.Abs(f(x[k + 1], y[k + 1]) - Math.Abs(f(x[k], y[k]))) < ee))
CurFlag = true;
}
if (CurFlag && PrevFlag)
break;
PrevFlag = CurFlag;
x[k] = x[k + 1];
y[k] = y[k + 1];
grad[0] = dfdx(x[k], y[k]);
grad[1] = dfdy(x[k], y[k]);
counter++;
}
Hisse[0, 0] = dfdx2(x[k], y[k]);
Hisse[0, 1] = dfdxdy(x[k], y[k]);
Hisse[1, 0] = dfdydx(x[k], y[k]);
Hisse[1, 1] = dfdy2(x[k], y[k]);
Det = Hisse[0, 0] * Hisse[1, 1] - Hisse[0, 1] * Hisse[1, 0];
if (Hisse[0, 0] > 0 && Det > 0)
TrueResult = true;
list3.Add(x[k + 1], y[k + 1]);
if (TrueResult)
{
if (Math.Abs(y[k+1]) < ee)
y[k+1] = 0;
if (Math.Abs(x[k+1]) < ee)
x[k+1] = 0;
label6.Text = "x=" + Convert.ToString(x[k + 1]);
label7.Text = "y=" + Convert.ToString(y[k + 1]);
}
else
MessageBox.Show("Точка минимума не найдена");
}
public void Marquardt(double x0, double y0, double ee)
{
double m = 100;
double[] grad = new double[2];
double[] d = new double[2];
double[] x = new double[2];
double[] y = new double[2];
double[,] Hisse = new double[2, 2];
double[,] HisseBuf = new double[2, 2];
double Det, temp;
x[0] = x0;
y[0] = y0;
list2.Add(x[0],y[0]);
int k = 0;
grad[0] = dfdx(x[k], y[k]);
grad[1] = dfdy(x[k], y[k]);
while (Math.Sqrt(Math.Pow(grad[0], 2) + Math.Pow(grad[1], 2)) > ee)
{
Hisse[0, 0] = dfdx2(x[k], y[k]);
Hisse[0, 1] = dfdxdy(x[k], y[k]);
Hisse[1, 0] = dfdydx(x[k], y[k]);
Hisse[1, 1] = dfdy2(x[k], y[k]);
while(true)
{
HisseBuf[0, 0] = Hisse[0, 0] + m;
HisseBuf[1, 1] = Hisse[1, 1] + m;
HisseBuf[0, 1] = Hisse[0, 1];
HisseBuf[1, 0] = Hisse[1, 0];
Det = HisseBuf[0, 0] * HisseBuf[1, 1] - HisseBuf[0, 1] * HisseBuf[1, 0];
temp = HisseBuf[0, 0];
HisseBuf[0, 0] = HisseBuf[1, 1]/Det;
HisseBuf[1, 1] = temp/Det;
HisseBuf[0, 1] = -HisseBuf[0, 1]/Det;
HisseBuf[1, 0] = -HisseBuf[1, 0]/Det;
d[0] = -(HisseBuf[0, 0] * grad[0] + HisseBuf[0, 1] * grad[1]);
d[1] = -(HisseBuf[1, 0] * grad[0] + HisseBuf[1, 1] * grad[1]);
x[k + 1] = x[k] + d[0];
y[k + 1] = y[k] + d[1];
if(f(x[k + 1], y[k + 1]) < f(x[k ], y[k ]))
{
m=m/2;
break;
}
else
m = 2 * m;
}
x[k] = x[k + 1];
y[k] = y[k + 1];
grad[0] = dfdx(x[k], y[k]);
grad[1] = dfdy(x[k], y[k]);
list2.Add(x[k + 1], y[k + 1]);
}
bool TrueResult=false;
if (Math.Abs(y[k + 1]) < ee)
y[k + 1] = 0;
if (Math.Abs(x[k + 1]) < ee)
x[k + 1] = 0;
Hisse[0, 0] = dfdx2(x[k], y[k]);
Hisse[0, 1] = dfdxdy(x[k], y[k]);
Hisse[1, 0] = dfdydx(x[k], y[k]);
Hisse[1, 1] = dfdy2(x[k], y[k]);
Det = Hisse[0, 0] * Hisse[1, 1] - Hisse[0, 1] * Hisse[1, 0];
if (Hisse[0, 0] > 0 && Det > 0)
TrueResult = true;
list3.Add(x[k + 1], y[k + 1]);
if (TrueResult)
{
label6.Text = "x=" + Convert.ToString(x[k + 1]);
label7.Text = "y=" + Convert.ToString(y[k + 1]);
}
else
MessageBox.Show("Точка минимума не найдена");
}
public void Graph()
{
for (double i = -1.9; i < 1.9; i += 0.005)
for (double j = 1.9; j > -1.9; j -= 0.005)
{
if (Math.Abs(f(i, j)) % 0.03 < 0.003)
list1.Add(i, j);
}
LineItem myCurve = myPane.AddCurve("", list1, Color.Green, SymbolType.Circle);
myCurve.Symbol.Fill.Color = Color.Green;
myCurve.Line.IsVisible = false;
myCurve.Symbol.Fill.Type = FillType.Solid;
myCurve.Symbol.Size = 1;
LineItem myCurve2 = myPane.AddCurve("", list2, Color.Red, SymbolType.Circle);
myCurve2.Symbol.Fill.Color = Color.Red;
myCurve2.Line.IsVisible = true;
myCurve2.Symbol.Fill.Type = FillType.Solid;
myCurve2.Symbol.Size = 5;
LineItem myCurve3 = myPane.AddCurve("", list3, Color.Blue, SymbolType.XCross);
myCurve3.Symbol.Fill.Color = Color.Blue;
myCurve3.Line.IsVisible = false;
myCurve3.Symbol.Fill.Type = FillType.Solid;
myCurve3.Symbol.Size = 10;
myPane.XAxis.Scale.Min = -2;
myPane.XAxis.Scale.Max = 2;
myPane.YAxis.Scale.Min = -2;
myPane.YAxis.Scale.Max = 2;
}
public void Newton(double x0, double y0,double ee)
{
double t=1;
double[] grad = new double[2];
double[] d = new double[2];
double[] x = new double[2];
double[] y = new double[2];
double [,] Hisse=new double[2,2];
double Det, temp;
int flag;
bool CurFlag=false, PrevFlag=false, TrueResult=false;
x[0] = x0;
y[0] = y0;
list2.Add(x[0], y[0]);
int k = 0;
int counter = 0;
grad[0] = dfdx(x[k], y[k]);
grad[1] = dfdy(x[k], y[k]);
while (Math.Sqrt(Math.Pow(grad[0], 2) + Math.Pow(grad[1], 2)) > ee)
{
Hisse[0, 0] = dfdx2(x[k], y[k]);
Hisse[0, 1] = dfdxdy(x[k], y[k]);
Hisse[1, 0] = dfdydx(x[k], y[k]);
Hisse[1, 1] = dfdy2(x[k], y[k]);
Det = Hisse[0, 0] * Hisse[1, 1] - Hisse[0, 1] * Hisse[1, 0];
temp = Hisse[0, 0];
Hisse[0, 0] = Hisse[1, 1]/Det;
Hisse[1, 1] = temp/Det;
Hisse[0, 1] = -Hisse[0, 1]/Det;
Hisse[1, 0] = -Hisse[1, 0]/Det;
if (Hisse[0, 0] > 0 && Det > 0)
{
d[0] = -(Hisse[0, 0] * grad[0] + Hisse[0, 1] * grad[1]);
d[1] = -(Hisse[1, 0] * grad[0] + Hisse[1, 1] * grad[1]);
t = 1;
flag = 0;
}
else
{
d[0] = -grad[0];
d[1] = -grad[1];
flag = 1;
t = 1;
}
x[k + 1] = x[k] + t * d[0];
y[k + 1] = y[k] + t * d[1];
if (flag == 1)
{
double boool = f(x[k], y[k]);
while (f(x[k + 1], y[k + 1]) >= boool)
{
t = t / 2;
x[k + 1] = x[k] + t * d[0];
y[k + 1] = y[k] + t * d[1];
}
}
list2.Add(x[k + 1], y[k + 1]);
if (counter == 0)
{
if ((Math.Sqrt(Math.Pow(x[k + 1] - x[k], 2) + Math.Pow(y[k + 1] - y[k], 2)) < ee) && (Math.Abs(f(x[k + 1], y[k + 1]) - Math.Abs(f(x[k], y[k]))) < ee))
PrevFlag = true;
}
else
{
if ((Math.Sqrt(Math.Pow(x[k + 1] - x[k], 2) + Math.Pow(y[k + 1] - y[k], 2)) < ee) && (Math.Abs(f(x[k + 1], y[k + 1]) - Math.Abs(f(x[k], y[k]))) < ee))
CurFlag = true;
}
if (CurFlag && PrevFlag)
break;
PrevFlag = CurFlag;
x[k] = x[k + 1];
y[k] = y[k + 1];
grad[0] = dfdx(x[k], y[k]);
grad[1] = dfdy(x[k], y[k]);
counter++;
}
Hisse[0, 0] = dfdx2(x[k], y[k]);
Hisse[0, 1] = dfdxdy(x[k], y[k]);
Hisse[1, 0] = dfdydx(x[k], y[k]);
Hisse[1, 1] = dfdy2(x[k], y[k]);
Det = Hisse[0, 0] * Hisse[1, 1] - Hisse[0, 1] * Hisse[1, 0];
if (Hisse[0, 0] > 0 && Det > 0)
TrueResult = true;
list3.Add(x[k + 1], y[k + 1]);
if (TrueResult)
{
if (Math.Abs(y[k+1]) < ee)
y[k+1] = 0;
if (Math.Abs(x[k+1]) < ee)
x[k+1] = 0;
label6.Text = "x=" + Convert.ToString(x[k + 1]);
label7.Text = "y=" + Convert.ToString(y[k + 1]);
}
else
MessageBox.Show("Точка минимума не найдена");
}
2
Документ
Категория
Рефераты
Просмотров
58
Размер файла
577 Кб
Теги
курсовик, отчет
1/--страниц
Пожаловаться на содержимое документа