Исследование методов решения линейных уравнений

Метод ветвей и границ

Метод первого порядка

Метод второго порядка

Оптимальный поиск

Пассивный поиск

метод дихотомии

Метод золотого сечения

Динамическое программирование

Линейное программирование. Метод Жордана-Гаусса.

Метод ветвей и границ

В качестве примера метода ветвей и границ рассмотрим задачу коммивояжера.

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

Имеется N пунктов, расстояние между которыми известны. Необходимо, начав движения из п. 1. последовательно обойти все остальные по самому короткому пути и вернуться в п. 1. Условие заключается в следующем: в каждый пункт надо войти только один раз и один раз из него выйти. На рис. 1. приведен пример задачи для N = 5. На ребрах, соединяющих вершины графа (пункты), указаны расстояния.

Рис. 1.

Условие делает невозможным выбор на некотором шаге вычислений оптимального из путей, приводящих в некоторую точку. Решение задачи методом прямого перебора всех возможных (N-1)! возможно лишь для малых N, Необходим метод, который позволяет уменьшать число рассматриваемых вариантов. Таким методом является метод ветвей и границ (МВТ).

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

Решение задачи.

Запишем исходные данные в виде таблицы C0 = (Cij). В ней i и j — номера пунктов, движение на каждом шаге совершается из i в j, длина этого шага пути Cij, крестиками отмечены запрещенные переходы. Пусть Сi = min Cij. Найдем Сi для каждого i (в таблице эти значения подчеркнуты), после чего совершим приведение матрицы С0 к C01 по строкам, где Cij1 вычисляются по формулам Cij1 = Cij — Ci. Обозначим d01 = Ci. и Cj1 = min Cij. Найдем Cj для всех j и выполним приведение C01 по столбцам. В приведенном примере все Сj = 0, поэтому C02 = C01, константа приведения по столбцам d02 = 0, d0 = d01 + d02 = 23.

Можно утверждать, что задачам с матрицами C0 и C01 соответствует один и тот же оптимальный маршрут. Длина любого пути LS (C0) и LS (C02) связаны формулой LS (C0) = LS (C02) + d0. что следует из условия. Тогда любой путь по матрице представляет движение по клеткам таблицы, причем в каждой строке и каждом столбце можно побывать только один раз.

Величину d0 можно использовать в качестве нижней границы при оценке всех возможных путей.

Рассмотрим путь, который включает переход из i в j. Для него LS ij = d0 + Cij. Такой переход всегда есть, так как в приведенной матрице хотя бы один элемент в строке — нулевой. Выбирая один из них, мы выбираем оптимизацию одного шага пути (самый короткий первый шаг). При этом, конечно, мы не знаем, как это отразится на длине последующего пути. Обозначим (ij) множество путей, не содержащих непосредственный переход из i в j. Поскольку из i надо куда-то выйти, а в j откуда-то прийти, то этому множеству соответствует оценка:

LS (ij) = d0 + ij, где ij =

Тогда нижней оценкой для множества путей будет 2 = d0 + ij. Мы заинтересованы в выборе такого перехода ij, для которого эта оценка является самой высокой. Этот выбор можно сделать, отыскав среди нулевых элементов i строки такой, которому соответствует 2 = max ij. Выбором пары (ij) все множество возможных путей разбивается на два подмножества. Одно из них содержит переход (ij) (ему соответствует оценка 1 = d0 + 1), другое (ij) (ему соответствует оценка 2 = d0 + 2).

В нашем примере первым шагом пути должен быть переход из п. 1. Сij= 0 для i, j = 1,2. Множество всех возможных путей разбиваем на два: содержащих переход (12) и не содержащих переход (12). Строим для них нижние оценки. 1 = d0 = 23. Для перехода, не содержащего (12), соответствует путь, содержащий переходы (14) и (32). На этом первая итерация вычислительного процесса заканчивается. При дальнейшем движении за поиском оптимального пути будем наблюдать с помощью дерева вариантов (рис.2).

Рис. 2.

Рассмотрим теперь первое подмножество. В нем уже реализован переход (12), поэтому в таблице С02 строку 1 и столбец 2 удаляем (вычеркиваем) из дальнейшего рассмотрения. С02 преобразуется в С1, а после приведения — в С12. В верхнем левом углу таблицы С1 ставим крестик, так как переход (21) в рассматриваемом подмножестве (12) становиться невозможным.

Теперь над новой таблицей проделываем те же действия, что и над С2. Для С12 константа приведения 1 = 2. Далее в качестве возможных рассмотрим подмножества, содержащие пути (23) и (), и произведем их оценку:

(23) 23 = 1 = (d0 + d1) + с23 = 23 + 2 + 0 = 25.

() = 2 = (d0 + d1) + 23 = 25 + 5 = 30,

где 23 = = C24 + C 43 = 3+2 = 5.

После вычеркивания из таблицы C12 строки i = 2 и столбца j = 3 получается матрица С2 размером 3×3. Очередное разбиение возможных вариантов перемещения на подмножестве (34) и () дает оценки 34 = 26 и =27.

После очередного сокращения матрицы таблица принимает размер 2×2 и однозначно определяет путь: начальная точка 4, конечная — 1, путь (451). Он имеет длину 6. По дереву вариантов мы прошли путь до вершины. Длинный путь (123 451) имеет длину 32. Аналогично множество () содержит единственный путь (3541) длиной 2. Путь (123 541) имеет длину 27, его длина меньше, чем у пути, содержащего переход (34). В качестве нижней границы теперь следует выбирать 27. Любое множество, имеющее оценку больше 27, следует отбросить. Если же оценка меньше 27, то такой путь подлежит исследованию. Если, пройдя по новому пути, получим значение длины меньше 27, то это новое значение следует взять в качестве новой оценки. Итак, мы прошли по дереву вариантов до вершины и длину одного их путей используем для последующей оценки вариантов.

Возвращаемся теперь по дереву вверх до ближайшего узла (12). Не рассмотренная ветвь дерева () имеет нижнюю оценку 30, что хуже достигнутого результата. Значит, эту ветвь можно далее не исследовать. Переходим к ветви (). Оценка для этой ветви ниже 27, поэтому данное множество следует рассмотреть.

В результате исследования придем к результату, что оптимальный путь (123 541).

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

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

Листинг программы

var

Way:TVector;

LengthOfWay:real;

HightLimitWasSeted:boolean;

function Addiction (var Matrix: Tarray;SizeOfMatrix:integer):real;

var

MinInRow, MinInCol: real;

i, j: integer;

d1,d2,d:real;

function MinElInRow (NumRow:integer):real;

var

i:integer;

min:real;

index:integer;

begin

index:=0;min:=1e38;

for i:=1 to SizeOfMatrix do

if (Matrix[NumRow, i]

min:=Matrix[NumRow, i]; index:=i

end;

if index=0 then raise Error;

MinElInRow:=min

end;

function MinElInCol (NumCol:integer):real;

var

i:integer;

min:real;

index:integer;

begin

index:=0;min:=1e38;

for i:=1 to SizeOfMatrix do

if (Matrix[i, NumCol]

min:=Matrix[i, NumCol]; index:=i

end;

if index=0 then raise Error;

MinElInCol:=min

end;

begin

d1:=0;

for i:=1 to SizeOfMatrix do

begin

MinInRow:=MinElInRow (i);

d1:=d1+MinInRow;

for j:=1 to SizeOfMatrix do

Matrix[i, j]: =Matrix[i, j]-MinInRow;

end;

d2:=0;

for i:=1 to SizeOfMatrix do

begin

MinInCol:=MinElInCol (i);

d2:=d2+MinInCol;

for j:=1 to SizeOfMatrix do

Matrix[j, i]: =Matrix[j, i]-MinInCol;

end;

d:=d1+d2;

Addiction:=d

end;

procedure FindBestWay (var matrix: TArray;SizeOfMatrix:integer;basis:integer;LowLimit:real);

var

LessMatrix:Tarray;

LowLimit1,LowLimit2:real;

i:integer;

d:real;

NextPoint:integer;

SpareDistance:real;

Row, Column: integer;

min1,min2,min:extended;

procedure Print;

var

i, j: integer;

number:integer;

s:string;

begin

for i:=1 to SizeOfMatrix do

begin

for j:=1 to SizeOfMatrix do

begin

if Matrix[i, j]>Maxint then s:=' X '

else Str (Round (Matrix[i, j]):3,s);

Text:=Text+s+' '

end;

Text:=Text+#13+#10

end;

Text:=Text+#13+#10

end;

procedure FindNextPoint (var BestPointAfterBasis: integer;var SpareDistance: real);

var

min1,min2:real;

i:integer;

RowOfBasis, BestColumn: integer;

begin

for i:=1 to SizeOfMatrix do

if Round (Matrix[i, 0])=basis then

begin

RowOfBasis:=i;

break

end;

for i:=1 to SizeOfMatrix do

if Matrix[RowOfBasis, i]=0 then

begin

BestColumn:=i;

BestPointAfterBasis:=Round (Matrix[0,i]);

break

end;

min1:=1e38;

for i:=1 to SizeOfMatrix do

if (i<>BestColumn) and (Matrix[RowOfBasis, i]

min1:=Matrix[RowOfBasis, i];

min2:=1e38;

for i:=1 to SizeOfMatrix do

if (i<>RowOfBasis) and (Matrix[i, BestColumn]

min2:=Matrix[i, BestColumn];

SpareDistance:=min1+min2;

end;

procedure FormLessMatrix (From, T0: integer);

var

i, j: integer;

ii, jj: integer;

RowOfFrom, ColumnOfTo, Row: integer;

begin

for i:=1 to SizeOfMatrix do

if Round (Matrix[i, 0])=From then

begin

RowOfFrom:=i;

break

end;

for i:=1 to SizeOfMatrix do

if Round (Matrix[0,i])=T0 then

begin

ColumnOfTo:=i;

break

end;

for i:=0 to SizeOfMatrix-1 do

for j:=0 to SizeOfMatrix-1 do

begin

if i>=RowOfFrom then ii:=i+1

else ii:=i;

if j>=ColumnOfTo then jj:=j+1

else jj:=j;

LessMatrix[i, j]: =Matrix[ii, jj]

end;

for i:=1 to SizeOfMatrix-1 do

if Round (LessMatrix[i, 0])=T0 then

begin

Row:=i;

break

end;

if Row<=SizeOfMatrix-1 then LessMatrix[Row, 1]: =1e38

end;

begin

Print;

if SizeOfMatrix=2 then

begin

min1:=Matrix[1,2]+Matrix[2,1];

min2:=Matrix[1,1]+Matrix[2,2];

if min1

else min:=min2;

LengthOfWay:=LowLimit+min;

if LengthOfBestway>LengthOfWay then

begin

inc (Way.CountOfPoints);

Way.Points[Way.CountOfPoints]: =Round (Matrix[0,2]);

LengthOfBestWay:=LengthOfWay;

BestWay:=Way;

dec (Way.CountOfPoints);

end;

if not (HightLimitWasSeted) then

begin

HightLimitWasSeted:=true;

HightLimit:=LengthOfWay

end;

end

else if SizeOfMatrix>2 then

begin

inc (Way.CountOfPoints);

repeat

FindNextPoint (NextPoint, SpareDistance);

Way.Points[Way.CountOfPoints]: =NextPoint;

FormLessMatrix (Basis, NextPoint);

d:=Addiction (LessMatrix, SizeOfMatrix-1);

LowLimit1:=d+LowLimit;

LowLimit2:=LowLimit+SpareDistance;

FindBestway (LessMatrix, SizeOfMatrix-1,NextPoint, LowLimit1);

if LengthOfBestWay>=LowLimit2 then

begin

for i:=1 to SizeOfMatrix do

if Round (Matrix[i, 0])=Basis then

begin

Row:=i;

break

end;

for i:=1 to SizeOfMatrix do

if Round (Matrix[0,i])=NextPoint then

begin

Column:=i;

break

end;

Matrix[Row, Column]: =1e38;

Addiction (Matrix, SizeOfMatrix);

LowLimit:=LowLimit2

end

until LengthOfBestWay<=LowLimit2;

dec (Way.CountOfPoints);

end

else raise Error

end;

begin

LengthOfBestWay:=1e38;

Way.CountOfPoints:=1;

Way.Points[1]: =1;

Text:='';

d0:=Addiction (MatrixOfLinking, n);

HightLimitWasSeted:=false;

FindBestWay (MatrixOfLinking, n,1,d0)

end;

end.

Метод первого порядка

В этом методе Xk+1=Xk-hkf (Xk). Очередная точка выбирается в направлении скорейшего убывания функции. Шаг hk должен выбираться достаточно малым, чтоб выполнялось f (Xk+1)< f (Xk)(рис 3).

Рис. 3 Рис.4

Величина шага должна уменьшаться по мере приближения к X0. тогда у нас происходит дробление шага. Шаг hk выбирают так, чтоб выполнялось неравенство

f (Xk+1) — f (Xk)? — е hk// f (Xk)//2

где 0 < е < 1 — произвольно выбранная константа.

На очередном шаге процесса проверяется неравенство. Если оно выполнено — шаг сохраняется, если нет — шаг уменьшается так, чтобы выполнилось неравенство. Геометрическая интерпретация представлена на рис. 4. траектория движения аппроксимирует ломаной линией градиентную кривую, проходящую через точки X0 и X0. Если выбрана постоянная величина шага, то количество итераций, которое необходимо выполнить, может оказаться очень большим.

Разновидностью является метод покоординатного спуска (рис.5). Движение к X0 осуществляется по отрезкам прямых, параллельным координатным осям.

Градиентные методы очень просты в реализации и поэтому часто используются на практике. Но при сложной структуре f (Xk) они могут плохо сходится. Поэтому используют вторые производные f (Xk).

Рис.5

Метод второго порядка

Рис.6

Методы второго порядка являются обобщенным методом Ньютона. На рис. 6 приведена геометрическая интерпретация отыскания корня уравнения (X)=0 методом касательных. Разложение (X) в окрестности Xk дает

(X) = (Xk) + (Xk)(X — Xk) + 0(| X — Xk |)

Отбрасывая малые высшего порядка и пологая, что (X)=0,получим

X = Xk+1 = Xk — (Xk) / (Xk).

Пусть теперь (X) — функция n-мерного аргумента. Разложим ее в ряд Тейлора, отбросив малые высшего порядка

(Xk) + (Xk)(X — Xk)., откуда Xk+1 = Xk — (Xk) / (Xk).

Теперь считаем (X)градиентом f (X), т. е. = f. Приравнивая f=0, находим стационарные точки f (X). Формула для вычисления координат этих точек преобразуется к виду:

Xk+1 = Xk — f (Xk) / f (Xk).

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

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

Листинг программы по методам 1го и 2го порядков

function f (x:real):real;

begin

// result:=sin (x);

result:=(x-3)*(x-3)*(x-3)*(x-3);

end;

function df (x:real):real;

begin

// result:=cos (x);

result:=4*(x-3)*(x-3)*(x-3);

end;

function ddf (x:real):real;

begin

//result:=-sin (x);

result:=12*(x-3)*(x-3);

end;

procedure TForm1. Button1Click (Sender: TObject);

var e, h, xk, xk1, d:real;

i:integer;

begin

e:=StrToFloat (Edit3.Text);

xk1:=StrToFloat (Edit1.Text);

h:=StrToFloat (Edit4.Text);

Memo1.Clear;

Memo1.Text:='Поиск с использованием f''(x)'#13#10#13#10;

i:=0;

repeat

xk:=xk1;

repeat

xk1:=xk-h*df (xk);

if df (xk)*df (xk1)>0 then h:=h*2;

if f (xk)

until (df (xk)*df (xk1)<=0) and (abs (f (xk))>=abs (f (xk1)));

Memo1.Lines.Add (Format ('X (%d)= %4.4f f (%1:4.4f) = %2:4.4f f''(%1:4.4f)=%3:4.4f h=%4:4.4f',[i, xk, f (xk), df (xk), h]));

Inc (i);

until (abs (xk1-xk)

Memo1.Lines.Add (Format (#13#10'Найдено значение f (%4.4f)=%4.4f с погрешностью e=%4.4f',[xk, f (xk), e]));

end;

procedure TForm1. Button2Click (Sender: TObject);

var e, h, xk, xk1, d:real;

i:integer;

begin

e:=StrToFloat (Edit3.Text);

xk1:=StrToFloat (Edit1.Text);

h:=StrToFloat (Edit4.Text);

Memo1.Clear;

Memo1.Text:='Поиск с использованием f''(x)'#13#10#13#10;

i:=0;

repeat

xk:=xk1;

repeat

xk1:=xk-h*df (xk)/ddf (xk);

if df (xk)*df (xk1)>0 then h:=h*2;

if f (xk)

until (df (xk)*df (xk1)<=0) and (f (xk)>=f (xk1));

Memo1.Lines.Add (Format ('X (%d)= %4.4f f (%1:4.4f) = %2:4.4f f''(%1:4.4f)=%3:4.4f h=%4:4.4f',[i, xk, f (xk), df (xk), h]));

Inc (i);

until (abs (xk1-xk)

Memo1.Lines.Add (Format (#13#10'Найдено значение f (%4.4f)=%4.4f с погрешностью e=%4.4f',[xk, f (xk), e]));

end;

end.

Оптимальный поиск

функция программирование задача

Рассмотрим класс функций одного аргумента, которые называются унимодальными. Определим унимодальную функцию y = f (X) следующим образом:

1. f (X) задана на отрезке [a, b]. 2. пусть X1 и X2 е [a, b], причем X1 < X2 пусть X0 — точка, в которой f (X) достигает минимума (максимума); тогда из X1 > X0 следует f (X2) > f (X1), а из X2 < X0 следует f (X1) > f (X2).

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

рис. 7

минимум f (X) может находится во внутренней точке [a, b] или на границе. Первоначально о положении X0 ничего не известно, кроме того, что, X0 е [a, b]. [a, b] можно назвать отрезком неопределенности.

Выбор Xi и вычисление F (Xi) называется экспериментом. Во всех случаях после выполнения заданных экспериментов отрезок неопределенности становится уже. Последовательное повторение экспериментов позволит достигнуть величины отрезка неопределенности меньшей, чем любое наперед заданное е > 0. правило выбора Xi называется стратегией. Оптимальной называется стратегия, которая приводит к X0 (с точностью до е) за минимальное число шагов (экспериментов).

Обозначим длину отрезка неопределенности Ln (X1, X2,.. ., Xn), где n — число экспериментов;

f (Xk) = min f (Xi), 1? i? n,

где K — номер эксперимента, которому соответствует минимальное значение f (Xi). Тогда получим:

X2 — X1, k=1;

Ln (X1, X2,.. ., Xn, к) = Xk+1 — Xk-1, 1? k? n;

Xn — Xn-1, k=n.

Т.о., Ln зависит как от выбора способа разбиения [a, b], так и от поведения f (X).

Пассивный поиск

В пассивном поиске все эксперименты проводятся одновременно. Без ограничения общности положим X1 = a = 0; Xn = b = 1;

Рассмотрим n = 4; 0 < X2 < X3 < 1

Тогда Ln = max { X2 — 0, X3 — 0, 1 — X2, 1 — X3}= max { X3, 1 — X2}

1? k? 4

Ln = min max { X3, 1 — X2}

0 < X2 < X3 < 1 k= 2,3

Так как X2 < X3, то оптимальной стратегией будет X2 = Ѕ - Д /2, X3 = Ѕ + Д /2, а соответствующее Ln = Ѕ + Д /2 (рис. 8,а).

Под Д > 0 понимают минимальное число, которое делает 2 эксперимента X2 и X3 различными.

При рассмотрении n = 5 получаем (рис 8, б), что Ln улучшилось на Д /2, что заставляет отказаться от нечетного числа экспериментов.

При произвольном четном n наилучшая стратегия составляется парами Д — экспериментов (Xi разнесены на расстояние Д).

В общем случае Ln? 2/n + Д /2.

рис.8

метод дихотомии

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

задать начальный интервал [a;b];

убедиться, что на концах функция имеет разный знак

< 0.;

Начальное приближение

x0 = .

повторять выбрать внутри интервала точку X;

сравнить знак функции в точке X со знаком функции в одном из концов;

если совпадает, то переместить этот конец интервала в точку X,

иначе переместить в точку X другой конец интервала;

После каждой итерации отрезок, содержащий корень уменьшается вдвое;

пока не будет достигнута нужная точность.

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

Метод выбора точки деления — ключевой для скорости работы метода. Хорошо, если функция, чей корень находится, проста и быстро вычисляется; но внутри функции могут быть и матричные операции, и численное интегрирование, и все что угодно; так что главным критерием оптимизации является минимизация числа делений (естественно, без ухудшения точности метода).

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

Метод золотого сечения

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

Воспользуемся соотношением

Lj-1 = Lj + Lj+1.

Будем поддерживать постоянными отношения длин интервалов

Lj-1 / Lj = Lj / Lj+1 = ф.

Это уравнение имеет единственный положительный корень ф = (1+5¹/2) / 2= 1.618 (рис.9). первые 2 эксперимента Рис.9

находятся на расстоянии 0.618 от концов отрезка. По результатам экспериментов сохраняется 1 из интервалов, и все повторяется. После n экспериментов

Ln = 1/ фn-1.

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

Листинг программы по методу пассивного поиска, методу дихотомии, методу «золотого сечения», оптимальный поиску

function TMyForm. F (x:real):real;

begin

Result:=(2*x-2)*(2*x-2)+2

end;

function TMyForm. PasPoisk (XMax, XMin, Eps, D: Real):string;

var

N: word;

i: word;

Index: word;

temp: real;

x: array[word] of real;

begin

n:=trunc ((2*(XMax-XMin))/Epsilon);

if (frac ((2*(XMax-XMin))/Epsilon)<>0) then inc (N);

if odd (N) then inc (N);

x[1]: =XMin;

x[N]:=XMax;

temp:=XMin;

for i:=1 to (N div 2 — 1) do

begin

temp:=temp+(XMax-XMin)/(N/2);

x[2*i]: =temp-D;

x[2*i+1]:=temp+D

end;

Index:=1;

for i:=2 to N do

if (f (x[Index])>f (x[i])) then Index:=i;

for i:=1 to N do chGraph. Series[0]. AddXY (x[i], f (x[i]),'', clBlue);

chGraph.Series[1].AddXY (x[Index], f (x[Index]),'', clRed);

Result:=FloatToStr (x[Index])+'; Количество итераций = '+IntToStr (N)+'.'

end;

function TMyForm. Dihotomy (XMax, XMin, Eps, D: Real):string;

var

N: word;

x: array [1.4] of real;

begin

x[1]: =XMin;

x[4]:=XMax;

N:=2;

chGraph.Series[0].AddXY (x[1], f (x[1]),'', clAqua);

chGraph.Series[0].AddXY (x[4], f (x[4]),'', clAqua);

while ((abs (x[4]-x[1])/2)>Eps) do

begin

x[2]: =((x[4]-x[1])/2)-D;

x[3]:=((x[4]-x[1])/2)+D;

inc (N, 2);

chGraph.Series[0].AddXY (x[2], f (x[2]),'', clAqua);

chGraph.Series[0].AddXY (x[2], f (x[3]),'', clAqua);

if (f (x[3])>f (x[2])) then x[4]: =x[3]

else x[1]: =x[2]

end;

chGraph.Series[1].AddXY ((x[1]+x[4])/2,f ((x[1]+x[4])/2),'', clRed);

Result:=FloatToStr ((x[1]+x[4])/2)+'; Количество итераций = '+IntToStr (N)

end;

function TMyForm. GoldSection (XMax, XMin, Eps, D: Real):string;

const Tau=0.618;

var

i: byte;

N: word;

x: array[1.4] of Real;

y: array[1.4] of Real;

begin

x[1]: =XMin;

x[4]:=XMax;

x[2]:=XMax-Tau*(XMax-XMin);

x[3]:=XMin+Tau*(XMax-XMin);

for i:=1 to 4 do y[i]: =f (x[i]);

for i:=1 to 4 do chGraph. Series[0]. AddXY (x[i], f (x[i]),'', clRed);

N:=4;

while (abs (x[4]-x[1])/2)>Eps do

if y[3]>y[2] then

begin

x[4]: =x[3]; y[4]: =y[3];

x[3]:=x[2]; y[3]: =y[2];

x[2]:=x[4]-Tau*(x[4]-x[1]);

y[2]:=f (x[2]);

chGraph.Series[0].AddXY (x[2], y[2],'', clRed);

inc (N)

end

else

begin

x[1]:=x[2]; y[1]: =y[2];

x[2]:=x[3]; y[2]: =y[3];

x[3]:=x[1]-Tau*(x[4]-x[1]);

y[3]:=f (x[3]);

chGraph.Series[0].AddXY (x[3], y[3],'', clRed);

inc (N)

end;

chGraph.Series[1].AddXY ((x[1]+x[4])/2,f ((x[1]+x[4])/2),'', clBlue);

Result:=FloatToStr ((x[4]+x[1])/2)+'; Количество итераций = '+IntToStr (N)

end;

function TMyForm. OptPoisk (XMax, XMin, Eps, D: Real):string;

var fb: array[byte] of word;

x: array[1.4] of real;

L: array[byte] of real;

i: byte;

N: word;

begin

x[1]: =XMin;

fb[0]:=1;

x[4]:=XMax;

fb[1]:=1;

N:=1;

while (fb[N]<((XMax-XMin)/Eps)) do

begin

inc (N);

fb[N]: =fb[N-1]+fb[N-2]

end;

L[N]:=((XMax-XMin)-fb[N-2]*D)/fb[N];

for i:=1 to N-1 do

L[i]: =fb[N-i+1]*L[N]-fb[N-i-1]*D;

for i:=1 to N do

begin

x[2]: =x[4]-L[i];

x[3]:=x[1]+L[i];

chGraph.Series[0].AddXY (x[2], f (x[2]),'', clBlue);

chGraph.Series[0].AddXY (x[3], f (x[3]),'', clBlue);

if f (x[3])>f (x[2]) then

begin

x[4]: =x[3];

x[3]:=x[2]

end

else

begin

x[1]:=x[2];

x[2]:=x[3]

end;

end;

chGraph.Series[1].AddXY ((x[1]+x[4])/2,f ((x[1]+x[4])/2),'', clRed);

Result:=FloatToStr ((x[4]+x[1])/2)+'; Количество итераций = '+IntToStr (N)

end;

procedure TMyForm. ExitClick (Sender: TObject);

begin

close

end;

procedure TMyForm. btExecClick (Sender: TObject);

begin

try

XMin:=StrToFloat (edMin.Text);

XMax:=StrToFloat (edMax.Text);

Epsilon:=StrToFloat (edEps.Text);

delta:=0.25*Epsilon

except

ShowMessage ('Данные некорректны');

end;

edMin.Enabled:=False;

edMax.Enabled:=False;

edEps.Enabled:=False;

btExec.Enabled:=False;

rbM1.Enabled:=False;

rbM2.Enabled:=False;

rbM3.Enabled:=False;

rbM4.Enabled:=False;

MyForm.Width:=500;

MyForm.Height:=310;

chGraph.Show;

pnAnswer.Show;

if rbM1. Checked=true then lbX. Caption:=lbX.Caption+(PasPoisk (XMax, XMin, Epsilon, Delta));

if rbM2. Checked=true then lbX. Caption:=lbX.Caption+(Dihotomy (XMax, XMin, Epsilon, Delta));

if rbM3. Checked=true then lbX. Caption:=lbX.Caption+(GoldSection (XMax, XMin, Epsilon, Delta));

if rbM4. Checked=true then lbX. Caption:=lbX.Caption+(OptPoisk (XMax, XMin, Epsilon, Delta));

end;

end.

Динамическое программирование

Метод динамического программирования — широко известный и мощный математический метод современной теории управления.

Рассмотрим управляемую систему, состояние которой в каждый момент времени характеризуется n-мерным вектором x с компонентами x1, …, xn. Предполагаем, что время t изменяется дискретно и принимает целочисленные значения 0, 1, … Так, для процессов в экономике и экологии дискретным значениям времени могут отвечать дни, месяцы или годы, а для процессов в электронных устройствах интервалы между соседними дискретными моментами времени равны времени срабатывания устройства. Предполагаем, что на каждом шаге на систему оказывается управляющее воздействие при помощи m-мерного вектора управления u с компонентами u1, …, um. Таким образом, в каждый момент времени t состояние системы характеризуется вектором x (t), а управляющее воздействие — вектором u (t). На выбор управления обычно бывают наложены ограничения, которые в достаточно общей форме можно представить в виде

u (t) k U, t = 0, 1, …

Здесь U — заданное множество в n-мерном пространстве.

Под влиянием выбранного в момент t управления (принятого решения) система переходит в следующий момент времени в новое состояние. Этот переход можно описать соотношением

x (t + 1) = f (x (t), u (t)), t = 0, 1, …

Здесь f (x, u) — n-мерная функция от n-мерного вектора x и m-мерного вектора u, характеризующая динамику рассматриваемой системы. Эта функция предполагается известной (заданной) и отвечает принятой математической модели рассматриваемого управляемого процесса.

Зададим еще начальное состояние системы

x (0) = x0,

где x0 — заданный n-мерный вектор. Таким образом, многошаговый процесс управления описывается соотношениями (1)-(3). Процедура расчета конкретного процесса сводится к следующему. Пусть в некоторый момент t состояние системы x (t) известно. Тогда для определения состояния x (t + 1) необходимо выполнить две операции:

1) выбрать допустимое управление u (t), удовлетворяющее условию (1);

2) определить состояние x (t + 1) в следующий момент времени согласно (2). Так как начальное состояние системы задано, то описанную процедуру можно последовательно выполнить для всех t = 0, 1, … Последовательность состояний x (0), x (1), … часто называется траекторией системы.

Заметим, что выбор управления на каждом шаге содержит значительный произвол. Этот произвол исчезает, если задать цель управления в виде требования минимизации (или максимизации) некоторого критерия оптимальности. Таким образом мы приходим к постановке задачи оптимального управления.

Задача оптимального управления Пусть задан некоторый критерий качества процесса управления (критерий оптимальности) вида Здесь R (x, u) и F (x) — заданные скалярные функции своих аргументов, N — момент окончания процесса, N > 0. При этом функция R может отражать расход средств или энергии на каждом шаге процесса, а функция F — характеризовать оценку конечного состояния системы или точность приведения в заданное состояние.

Задача оптимального управления формулируется как задача определения допустимых управлений u (0), u (1), …, u (N — 1), удовлетворяющих ограничениям (1), и соответствующей траектории, то есть последовательности x (0), x (1), …, x (N), которые в совокупности доставляют минимальное значение критерию (4) для процесса (2), (3).

Минимизация критерия (4) обычно отвечает выбору управления, обеспечивающего наименьшие затраты средств, ресурсов, энергии, наименьшее отклонение от заданной цели или заданной траектории процесса. Наряду с этим часто ставится также задача о максимизации критерия вида (4), например о максимизации дохода или объема производства. Однако нетрудно видеть, что максимизация критерия J эквивалентна минимизации критерия (- J). Поэтому простая замена знака у функций R и F в (4) приводит задачу о максимизации критерия к задаче о его минимизации. Далее всюду для определенности рассматриваем задачу о минимизации критерия (4).

Листинг программы

type

ElType = string[3];

TVector = array of ElType;

const FORBIDDEN: ElType = 'x';

var

Form1: TForm1;

D: Array of TVector; // Матрица весов переходов

P: Array of TVector; // Матрица узлов

implementation

{$R *.dfm}

procedure TForm1. Button2Click (Sender: TObject);

var

i, j: integer;

_Width, _Height: Integer;

Weight1, Weight2: integer;

X, Y: Integer;

S: string;

begin

_Width := 5;

_Height := 4;

SetLength (D, _Width);

for i := Low (D) to High (D) do

SetLength (D[i], 2 * _Height — 1);

D[0,0] := '4'; D[1,0] := '5'; D[2,0] := '3';

D[0,1] := '5'; D[1,1] := '4'; D[2,1] := '6';

D[0,2] := '4'; D[1,2] := '5'; D[2,2] := '5';

D[0,3] := '4'; D[1,3] := '5'; D[2,3] := '6';

D[0,4] := '6'; D[1,4] := '5'; D[2,4] := '5';

D[0,5] := '7'; D[1,5] := '5'; D[2,5] := '7';

D[0,6] := '6'; D[1,6] := '8'; D[2,6] := '7';

D[3,0] := '5'; D[4,0] := FORBIDDEN;

D[3,1] := '6'; D[4,1] := '7';

D[3,2] := '5'; D[4,2] := FORBIDDEN;

D[3,3] := '6'; D[4,3] := '8';

D[3,4] := '7'; D[4,4] := FORBIDDEN;

D[3,5] := '8'; D[4,5] := '9';

D[3,6] := '7'; D[4,6] := FORBIDDEN;

SetLength (P, _Width);

for i := Low (P) to High (P) do

SetLength (P[i], _Height);

for i := Low (P) to High (P) do

for j := Low (P[i]) to High (P[i]) do begin

if (i=0) and (j=0) then begin

P[i, j] := '0';

Continue;

end;

if (i=0) then begin

P[i, j] := IntToStr (StrToInt (P[i, j-1]) + StrToInt (D[i, 2*j-1]));

Continue;

end;

if (j=0) then begin

P[i, j] := IntToStr (StrToInt (P[i-1,j]) + StrToInt (D[i-1,2*j]));

Continue;

end;

Weight1 := StrToInt (P[i, j-1]) + StrToInt (D[i, 2*j-1]);

Weight2 := StrToInt (P[i-1,j]) + StrToInt (D[i-1,2*j]);

if Weight1

P[i, j] := IntToStr (Weight1);

D[i-1,2*j] := FORBIDDEN;

end

else

begin

P[i, j] := IntToStr (Weight2);

D[i, 2*j-1] := FORBIDDEN;

end;

end;

for i := Low (P) to High (P) do begin

S := '';

for j := Low (P[i]) to High (P[i]) do

S := S + P[i, j] + #9;

Memo1.Lines.Add (S);

end;

X := _Width;

Y := _Height;

S := '(' + IntToStr (X) + ';' + IntToStr (Y) + ')';

while (X<>1) and (Y<>1) do begin

if D[X-1−1,2*Y-1−1]<>FORBIDDEN then X := X — 1

else Y := Y — 1;

S := '(' + IntToStr (X) + ';' + IntToStr (Y) + ') -> ' + S;

end;

S := '(1;1) -> ' + S;

ShowMessage (S);

end;

end.

Линейное программирование. Метод Джордана-Гаусса

Методы линейного программирования — численные методы решения оптимизационных задач, cводящихся к формальным моделям линейного программирования.

Как известно, любая задача линейного программирования может быть приведена к канонической модели минимизации линейной целевой функции с линейными ограничениями типа равенств. Поскольку число переменных в задаче линейного программирования больше числа ограничений (n > m), то можно получить решение, приравняв нулю (n — m) переменных, называемых свободными. Оставшиеся m переменных, называемых базисными, можно легко определить из системы ограничений-равенств обычными методами линейной алгебры. Если решение существует, то оно называется базисным. Если базисное решение допустимо, то оно называется базисным допустимым. Геометрически, базисные допустимые решения соответствуют вершинам (крайним точкам) выпуклого многогранника, который ограничивает множество допустимых решений. Если задача линейного программирования имеет оптимальные решения, то по крайней мере одно из них является базисным.

Метод Жордана-Гаусса (Jordan) почти ничем не отличается от классического метода Гаусса, только матрица системы приводится не к треугольному, а к диагональному единичному виду. Если отвлечься от перестановки строк и столбцов для выбора наивыгоднейшего ведущего элемента, то процесс выглядит так. Пусть нам дана система n линейных уравнений с таким же количеством неизвестных. Представим коэффициенты этой системы и свободные члены в виде расширенной матрицы A размером n*(n+1). Элементы матрицы будем обозначать a[i, j]. Выполним операции:

для каждого i от 1 до n: .

| для каждого k от n+1 до 1 с шагом (-1): .

| | a[i, k] <- a[i, k]/a[i, i] .

| для каждого j от 1 до n, не равного i: .

| | для каждого k от n+1 до i с шагом (-1): .

| | | a[j, k] <- a[j, k] - a[j, i]*a[i, k] .

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

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

Листинг программы

function Divide (var V: Vector; D: real): boolean;

var

I: Integer;

begin

Result := true;

for I := Low (V) to High (V) do

if (V[I]=0) and (D=0) then V[I] := 0 else

if D<>0 then

V[I] := V[I] / D

else Result := false;

end;

procedure Multiply (var V: Vector; K: real);

var

I: Integer;

begin

for I := Low (V) to High (V) do

V[I] := V[I] * K;

end;

procedure Subtract (var V1, V2: Vector);

var

I: Integer;

begin

for I := Low (V1) to High (V1) do

V1[I] := V1[I] - V2[I];

end;

procedure TForm1. Button2Click (Sender: TObject);

var

M: Matrix;

i, j: integer;

k: real;

Report: TStringList;

begin

SetLength (M, StringGrid1. RowCount-1);

for I := Low (M) to High (M) do

SetLength (M[I], StringGrid1. ColCount-1);

for I := Low (M) to High (M) do

for J := Low (M[I]) to High (M[I]) do

M[I, J] := StrToInt (StringGrid1.Cells[J+1,I+1]);

Report := TStringList. Create;

Report.Add ('

');

Report.Add ('

Данная СЛАУ:');

MatrixToHTML (Report, M);

for I := Low (M) to High (M) do begin

Report.Add ('

1. Избавляемся от коэффициентов перед x_{' + IntToStr (I+1) + '}:');

   if M[i, i]<>1 then

   if not Divide (M[I], M[I, I]) then begin // Приведение текущей строки

   MessageBox (Handle, 'Деление на ноль: …', 'Внимание', MB_OK + MB_ICONWARNING);

   Report.Free;

   Exit;

   end;

   for j := Low (M) to High (M) do begin

   if I=j then Continue;

   k := 1;

   if M[I, I]<>M[J, I] then begin

   k := M[J, I];

   if not Divide (M[J], k) then begin

   MessageBox (Handle, 'Деление на ноль: …', 'Внимание', MB_OK + MB_ICONWARNING);

   Report.Free;

   Exit;

   end;

   end;

   Subtract (M[j], M[I]);

   if (k<>1) and (I>j) then Multiply (M[j], k);

   end;

   MatrixToHTML (Report, M);

   end;

   Report.SaveToFile (ExtractFilePath (Application.ExeName) + 'result.html');

   Report.Free;

   MessageBox (Handle, PChar ('Создан файл отчета: ' + ExtractFilePath (Application.ExeName) + 'result.html'),

   'Решение получено', MB_OK + MB_ICONINFORMATION);

   WebBrowser1.Navigate (ExtractFilePath (Application.ExeName) + 'result.html');

   end;

   procedure TForm1. Button1Click (Sender: TObject);

   var

   N: Integer;

   begin

   N := StrToInt (LabeledEdit1.Text);

   StringGrid1.RowCount := N+1;

   StringGrid1.ColCount := N + 2;

   end;

   end.

   Заключение

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

   В результате выполнения работы были реализованы следующие алгоритмы:

   — Метод ветвей и границ на примере задачи о коммивояжоре;

   — Поиски 1-ого и 2-ого порядка;

   — Линейное программирование на примере симплексного метода.

   1. Галкин А. А. методы оптимизации в примерах и задачах: Учеб. Пособие. — Владимир: ВПИ, 1989. — 96 с.

   2. Жирков В. Ф. моделирование производственных процессов с дискретным характером управления: Учеб. Пособие. — Владимир: НПИ, 1984. — 84 с.

   3. Черноусько Ф. Л., Баничук Н. В. Вариационные задачи механики и управления: Численные методы. М.: Наука, 1973. 238 с.

   4. Internet.