ИЗВЕСТИЯ
ТОМСКОГО ОРДЕНА ОКТЯБРЬСКОЙ РЕВОЛЮЦИИ И ОРДЕНА ТРУДОВОГО КРАСНОГО ЗНАМЕНИ ПОЛИТЕХНИЧЕСКОГО ИНСТИТУТА ИМ. С! М. КИРОВА
Том 277 1977
ПРИМЕНЕНИЕ СЛУЧАЙНОГО ПОИСКА
В ОДНОЙ ЗАДАЧЕ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ
В. В. ЗАХАРОВ, В. П. ИВАНЧЕНКОВ, В. А. ФИЛИНОВА
(Представлена научным семинаром лаборатории вычислительной техники и автоматизации НИИ ЯФ ЭА)
Для решения задач нелинейного программирования небольшой размерности и с небольшим числом ограничений кажется естественным привлекать методы статистического поиска,, обеспечивающие простые машинные алгоритмы. В настоящей статье предлагается применить для этого метод случайного поиска [1] в рамках одной схемы функций штрафа [2].
1. Постановка задачи
• *
Пусть F(x) и fj (я) (/£/, /={1, 2, т}) — выпуклые функции векторного аргумента х={хи хп). В работе рассмотрена задача: найти
. min[F(X)\/j(X)>0, j£j\. (1)
X
Следует подчеркнуть, что сходимость рассматриваемого варианта метода штрафных функций доказана в [3] для выпуклых функций F, fj. Тем не менее в этой же работе были рассмотрены примеры решения невыпуклых задач «хорошего поведения», т. е. таких, в которых кривизна поверхностей F(x) —const, ff(x)= const не была слишком большой. Метод функций штрафа был выбран в настоящей работе в сочетании со случайным поисковым движением по той причине, что он позволял надеяться на приближенное решение рассматриваемой далее практической задачи, в отношении выпуклости которой не было уверенности.
2. Описание метода
Решение задачи в [2], [3] сводится к поэтапной минимизации функции
т]
Ат(Х)=2 Л(Х), (2)
представляющей собой штраф за приближение к границе области определения R функции F(X).
\
В настоящей работе на (¿+1)-м этапе поиска min Ат(Х) (об-ласть R = {fj(X)>0, / £Л) функция (2) дополнялась ограничением:
fm+i(X) = F(X°*)-F(X)>0,
где Xok —точка траектории поиска на k-м этапе, в которой функция Ат+г (Я) =А (X) достигает наименьшего значения.
Модификация использованного метода штрафа состояла еще в том, что на k-м этапе не ставилась задача отыскания min Л (Я) с большой точностью. Более того, выход из k-ro этапа минимизации происходил тогда, когда разрешающая способность метода случайного поиска (шаг метода) уменьшалась до некоторой фиксированной величины т.
3. Алгоритм
\
Рассмотрим алгоритмйческую реализацию метода на некотором, k-м этапе. Основной задачей этапа является получение точки такой, что F(X°'k+l ) <lF(Xok). Спуск к этой точке осуществляется по поверхности штрафа Л(А") методом случайного поиска с шагом т. Не трудно видеть, что fm+\ (Х0/г )=0 и А(Хок)—-\-оо. Поэтому неравенство для дополнительного (т+1)-го ограничения имеет смысл взять в виде
fm+l(X)=F(X**)-F(X)+Ф>0,
где ф>0.
Теперь уже можно искать ттЛ(Х) методом случайного поиска. В соответствии с известным алгоритмом случайного поиска с пересчетом [1] принимаем точку Xk за начальную точку траектории поиска на k-м этапе: Xk=Xo . В кубической окрестности точки (ребро куба равно т) Я"* производим равномерное рассеяние точек до тех пор, пока не найдется хотя бы одна из них, в которой значение Л(Х) меньше, чем Л(Хо). Обозначим эту удачную точку через Х\ , вновь примем ее за центр куба с ребром т и вновь будем производить равномерно распределенные испытания в этом кубе и т. д.
Для параметра т введем адаптацию к условиям поиска: величину т будем уменьшать в q раз, если нарушено хотя бы одно из (т+1) ограничений; величину т будем увеличивать в q раз в противном случае. Начальное значение полагалось равным т. Если выполнялось т<;т, то
£-й этап считался оконченным и полагалось х1=ХкАгХ , причем, очевидно,
А(Х*=Х%)> А (Х\)> • • • >A(XLi)>A{XkL=Xk+ly
Можно предположить (и практические вычисления это подтверждают), что на заключительных этапах поиск будет происходить в узком овраге. Вследствие этого в кубе с центром Xi необходимо провести чрезвычайно много испытаний, прежде чем выполнится А (Xf ) > >A(Xi+1 ). Против такой ситуации введено следующее видоизменение вычислительного алгоритма. Если количество испытаний в окрестности Х\ превысило (зацикливание), то из точки X* производится вторично спуск по описанному алгоритму до тех пор, пока не выполнился
тгСг (т. е. естественное окончание этапа) либо вновь произойдет за-
I ' & /? & *
цикливание в точке . Пусть А (Х/3) <Л (А"/, ). В этом случае в качестве Хь выбирается точка
V £ у^
А/?
- А -^-,
где точка Х% есть точка пересечения границы области # с прямой Х^ Х^ в направлении убывания А(Х) вдоль этой прямой. На этом к-й этап также оканчивается.
Вся процедура решения задачи (1) продолжается до тех пор, пока не будет исчерпано заданное заранее количество вычислений функции Р(Х), равное N.
4. Блок-схема и ее описание
с
Более детально строение алгоритма можно уяснить из приводимой блок-схемы. Поясним ее.
Блок 1 осуществляет формирование команд программы с учетом числа переменных п и числа ограничений т.
Блок 2. Признак й определяет номер спуска из начальной точки Х\ ; если 1, то из точки XI спуск совершен уже дважды. Рабочий шаг т в начале к-то этапа поиска принимает значение, равное т. В блоке 2 вычисляются начальные значения Р(Х) и А(Х).
Блок 3. Точка Хк — центр поиска на 6-м этапе. Начальное значение Хк равно Хо .
Блок 4. Счетчик (т) необходим для подсчета числа последовательных точек, в которых А (Хч) >Л (Хкг), где Xа —случайная, равномерно распределенная точка в кубе с центром
с с
Блоки 5—6. производят вычисление ¡] (X ). Если условие /ДА не выполняется, то шаг т уменьшается в ¿7 раз (блок 17).
Блок 7. Счетчик (И) контролирует число вычислений ^(Я), если cт{F)=N, то происходит печать найденных /7т1п и ХтЫ . Поиск прекращается (блок 18).
Блоки 8—9. М — текущее минимальное значение Р(Х)9 X' —точка, соответствующая М.
Блок 10. Рабочий шаг т увеличивается в ц раз, если для точки X выполнилось т ограничений.
Блок 11 проверяет: Л(Х^)<Л(Х? ). Если да, то-^блок 26: А{Х?^):=А(Х*)9Хи : = , иначе сг (т):- = сг (т)+1.
Блоки 12—16. Если сг (т)=Л^> то шаг т уменьшается в ц раз, иначе (в блоке 14) проверяется: с1=0?
Если да, то Х\х :=Хк и ¿¿=1.
С шагом т=т/2 еще раз осуществляется спуск до получения точки
XI .
Блоки 19—22. Если йФ0, тю : =
Делается проверка А (А*) <Л (А^г), чтобы выбрать направление движения по прямой X?1 Х\%.
1 Формирование команд, сч (^)=0
1
2 т :=т ' Вычисление
Ф
► з
1
> 4 сч (-): =0
1
> 5 ---1
ф нет
I
Вычисление /ДХ5)
/нет
-> 17
г : =х¡д, ?<т ?
да
I да
Вычисление ^(Х-), | да сч(Л:-сч(/р +11сч(/Г)-А^? |
18
Печать /Гш!п Останов.
I
нет
Р(Х*)<М?
| да
M:^F(X'), ХМ
Ф
10
1-Я
I
11
А(Х")<А(Х11)?
да
| нет
Л12
нет\_
сч(т) 2 =сч(т:)Ч-1 сч(х)^э?
ф да
19
Ф
20
А(Х1)<А{Х1) ?
/ да
\ нет
21 X? —нач. м
движ.
22 X?—нач. *а
движ.
Движение по прямой 23 х£ X* с шагом т/б
I
24
ЛГ^Х^хЯ?
£1 -о
'нет
I да
25
у/г | у&
у/г __ Л/,_1~ ЛЯ —о-
нет^З
т : — т/^, Х<-С ?
Ф да
14
нет
¿=0 ?
26
I да
15
1
16 т : —т/2
Блок-схема алгоритма
5 Заказ 10892
Блоки 23—26 осуществляют движение по прямой Хг\ Х\% с шагом т¡0, до границы области К, после чего определяется начальная точка следующего этапа поиска ттА(Х):
у к | -у к
5. Описание программы
Программа составлена для ЭВМ БЭСМ-4. а) параметры, используемые в методе.
Таблица
№ Ячейки Обозначение Значение в программе Назначение параметра
0042 10~5 Малая положительная величина
0043 с 1 Верхняя граница шага
0044 X 10-2 Нижняя граница шага
0045 Я 2 Коэффициент изменения шага
Максимально допустимое количест-
0046 N 103 во-вычислений функции цели
0047 А Большое положительное число
Коэффициент уменьшения шага при
0050 а 10 движении по прямой
Количество допустимых нарушений
0051 Хэ 2 неравенства Л(х^) < А(х Ке)
Эти параметры программист может менять по своему усмотрению, б) Распределение памяти и обращение к программе метода. Программа метода расположения в ячейках с 0030 по 0330. С 0400 по 1277 ячейку занимает рабочее поле для промежуточных и конечных результатов. Программа использует константы ИС-2.
Обращение* к программе осуществляется командами:
16
Х+\
0100
0077
Х+1 | 00 |
т
где
—число переменных, т — число ограничений на оптимизируемые параметры, Р — ячейка начала счета }¿{X) и Р{Х), п, т, И задаются в восьмеричном виде.
Программы вычисления ¡¿(Х) и Р(Х) могут быть размещены в любом месте памяти, начиная с 1300 ячейки. Ячейка /М предназначена для возврата в программу метода и должна быть свободной.
Компоненты вектора X занимают ячейки 0601—0600+п. Здесь же задает программист начальную точку Х°.
, Рекомендуется вначале вычислить (X) и занести эти значения в ячейки 2701—2700+т, затем проверить Если хотя бы
одно из ограничений не выполнилось, необходимо передать управление
в ячейку 0076, иначе вычислять F(X). Значение F(X) необходимо записать в. 0600 ячейку. По окончании счета F(X) передать управление в ячейку F—Л.
Текущее минимальное значение функции и точка, ему соответствующая, хранятся в ячейках 0400—0400+лг. Рабочий шаг i хранится в 0060 ячейке.
6. Пример применения метода к производственной задаче
Задача состояла в определении параметров тракта транспортировки заряженных частиц одного ускорителя, оптимальных с точки зрения минимума стоимости конструкции в целом. Тракт состоял из 13 участков: 7 квадрупольных линз и 6 участков дрейфа частиц. Каждый участок однозначно определяется набором параметров в [4]
li, Р/, «i Ti (¿=1,13) .
Величины Р;, aJt j¿ вычисляются по следующим формулам. Для горизонтальной плоскости фокусировки линз:
Pi=То/ sin2e¿-2a0i—=sin 0¿ cos ef+^cos20f, У k¡
аг=70г—=-cos вг sin 9¡+a0;(cos2 9£—sin20É) + kt sin 0¿cos 0¿, V k¿
Tí=ToíCos2 0j+2 a0¿l/lTcos 0¿sin 0¿+M¿ eos2 9t,(i = 1.....7). Для участков дрейфа пучка частиц:
ai—ao¿ То А
- «./
i — I 0¿
Для вертикальной плоскости фокусировки линз:
k = 0. (1/2=1,...,6).
Р^ТоЛсЬ3 в _2а0(4=сЬ е,съ в,+р0,сЬ2в0 а, = -ТпгТ7ксЬ 0г сИ 0,;+аО1.(сЬ20г+сЬ20г)-рог1/^сЬ 0Г сЬ 0,,
V к1
Тг = ТогсЬ20г-2аог> ^сЬ0гсЬ0г+^^сЬ2 в4,(/=1.....7).
Здесь в формулах всюду в; —/¿у к^ . В этих обозначениях функция цели ^(Я) имеет вид:'
7 _
/=1
с-С=const.
здесь I — номер участка тракта с линзой, ргор.л, Рверт.г —максимальные на участке г в вертикальной или горизонтальной плоскости фокусировки линз. Переобозначим вектор переменных
Х= ЛГ2 = &4, —
С "Г
На переменные х были наложены ограничения:
0,09^x^0,64; 0,09^х2^1; 4^х3<10;
4<хА<10; 4<х5^9; 4^хб^9;
2*1=31,92.
¿=з
Остальные параметры не менялись.
Для действующего тракта значения переменных были равны:
XI = £2=0,438; х3=8; х4=7,96;
^5=7,96; х6=,8; ^(Х) =0,08230.
После минимизации Р(Х) изложенным выше методом получили:
. Р(Х) =0,00135;
= 0,184; х2=0,375; х3 = 5,084;
х4 = 8,460; х5=8,895; х6 = 9,480.
Машинное время, затраченное на поиск /•'(Х*) =0,00135 из данной начальной точки, составило примерно 13 мин. Подходящий результат методом случайного поиска без всяких модификаций ранее был получен за 4 часа работы БЭСМ-4.
Алгоритм метода в том виде, как он приведен на блок-схеме, довольно чувствителен к оврагам, поэтому когда к имевшимся ограничениям были добавлены еще ограничения ^0,0833; 77(^)^33,3, имевшие ярко выраженный овражный вид, то метод быстро зацикливал.
Таким образом, для повышения эффективности метода следует увеличить его разрешающую способность в овражной ситуации.
Приложение
Программа метода
0100 56 - 0303 0056 1 1 32 0001 0130 7777
1 4 55 - 7734 0072 2 00 — — 0062
2 1 14 0064 0072 0072 3 72 — 0072 -
3 4 55 - 7732 '0073 4 54 0107 0040 0053
4 4 55 - 7731 0074 5 07 0053 0040 0053
5 14 0114 0074 0074 6 13 0041 0053 0053
6 00 0047 - 0400 7 23 0053 — 0040
7 00 — - 0067 0140 05 7762 0040 0054
0110 00 - - 0061 1 02 0054 7761 0054
1 72 - 0072 — 2 05 0054 00,60 0054
2 5 00 0600 — 0500 3 3 01 0054 0700 0600
3 1 32 0002 0112 7777 4 1 32 0002 0134 7777
4 00 -. — 0063 5 72 — 0073 -
5 00 0043 — 0060 6 3. 16 0147 - 7777
6 72 - 0073 — 7 01 7761 0067 0067
7 3 16 0120 - 7777 0150 15 0067 0046 -
0120 56 0600 0300 0500 1 36 — 0272 -
1 2 04 7761 2700 0075 2 02 0600 0400 -
2 01 0075 0065 0065 3 76 — 0157 0070
3 1 32 0002 0121 7777 4 72 - 0072 -
4 04 7761 0042 0075 5 5 00 0600 — 0400
5 01 0075 0065 0065 6 1 32 0001 0155 7777
6 00 0065 — 0064 7 05 0060 0045 0060
7 72 — 0072ч — 0160 72 •— 0074 —
0130 5 00 0500 — 0700 1 2 04 7761 2700 0075
2 01 0075 0070 0070 0220 00 — — —
3 1 32 0002 0161 7777 1 72 — 0072 .—
4 02 0700 0600 0075 2 5 00 0700 — 1100
5 01 0075 0042 0075 3 1 32 0002 0222 7777
6 56 — 0312 — 4 04 0060 0050 0060
7 01 0075 0070 0070 5 02 0066 0064 —
0170 02 0070 0064 — 6 36 - 0235 —
1 76 — 0177 — 7 00 - - —
2 72 — 0072 — 0230 72 -- 0072 —
3 5 00 0600 — 0700 1 4 00 1000 - 0075
4 1 32 0001 0173 7777 2 5 00 1100 - 1000
5 00 0070 — 0064 3 1 00 0075 - 1100
6 56 — 0132 — 4 1 32 0002 0231 7777
7 01 0062 7761 0062 5 72 - 0072 -
0200 15 0062 0051 — 6 5 00 1000 -- 0600
1 36 — 0203 — 7 5 01 0600 0060 0600
2 56 — ■ 0133 — 0240 6 02 0600 1000 0075
3 04 0060 0045 0060 1 6 02 1077 0777 0057
4 02 0060 0044 — . 2 05 0057 0075 0057
5 36 —: 0207 — 3 6 02 1100 1000 0075
6 56 — 0132 —■ 4 04 0057 0075 0075
7 15 0063 — — 5 3 01 0075 0777 0577
0210 76 — 0221 — 6 1 32 0003 0240 7777
1 00 7761 — 0063 7 00 0310 -- 0076
2 72 — 0072 — 0250 ' 16 0251 0145 0157
3 5 00 0700 — 1000 1 00 _ _ _
4 1 32 0002 0213 7777 2 72 _ 0072 __
5 00 0064 — 0066 3 5 00 0600 - 1200
6 04 0043 7762 0060 4 1 32 0002 0253 7777
7 56 — 0126 — 5 72 - 0072 -
6 56 — 0237 — 3 36 - 0266 -
7 00 0311 — 0157 4 04 7761 0075 0075
0260 72 —■ 0072 — 5 01 7761 0056 0056
1 6 01 1000 1200 0075 6 '56 _ 0167 _
2 1 04 0075 7762 0600 7 15 1001 0701 _
3 1 32 0002 0261 7777 0320 36 _ 0323 _
4 01 0061 7761 0061 1 72 _ 0072 _
5 56 0307 Olli 0076 2 56 _ 0222 _
6 04 0060 0045 0060 3 04 0060 7760 0060
7 02 0060 0044 — 4 56 — 0126 —
0270 76 — 0132 —
1 56 — 0207 — '
2 72 — 0072 —
3 16 0274 7501 7610
4 1 72 0400 0027 0400
5 16 0276 7501 7610
6 72 0060 0027 0063
7 77 — - -
0300 72 — 0074 -
1 56 - 0121 0065
0302 00 - — —
3 16 0304 7501 7610
4 52 0042 0042 0051
5 72 - 0077 -
6 56 0307 0101 0076
7 16 - 0266 —
0310 16 - 0257 —
1 05 0060 0045 0060
2 02 0075 - —
ЛИТЕРАТУРА
1. Jl. А. Растр.игин. 'Статистические методы поиска. М., «Наука», 1968.
2. С. W. Carroll. The CRST for Optimizing Nonlinear Restraines systems. Oprns. Res., v. 9, № 2, 1961.
3. A. V. Fiacco, G. P. Mc С о r m i с k. Computational Algorithm for SUMT for Nonlinear programmung. Manag. Sc., v. 10, № 4, 1964.
4. В. П. Иванченков. Вопросы применения корреляционно-экстремальных систем в задачах электронной оптики. Диссертация. Томск, 1969.