Модель сферического маятника с подвижной точкой подвеса в задаче пространственного перемещения груза грузоподъемным краном при ограничении колебаний Текст научной статьи по специальности «Физика»

14. Wu X., He X., Sun N. An analytical trajectory planning method for underactuated overhead cranes with constraints // Proceedings of the 33rd Chinese Control Conference, Nanjing. 2014. P. 1966-1971. DOI: 10.1109/ChiCC.2014.6896931.

15. Korytov M., Shcherbakov V., Titenko V. Analytical solution of the problem of acceleration of cargo by a bridge crane with constant acceleration at elimination of swings of a cargo rope // Journal of Physics: Conference Series. 2018. Vol. 944, No. 1. P. 012062. DOI: 10.1088/1742-6596/944/1/012062.

16. Блехман И. И. Вибрационная механика. М. : Физматлит, 1994. 400 с.

17. Shustov V. V. Approximation of functions by asymmetric two-point Hermite polynomials and its optimization // Computational mathematics and mathematical physics. 2015. Vol. 55, no. 12. P. 1960-1974.

18. Shustov V. V. Approximation of functions by two-point Hermite interpolating polynomials // Computational mathematics and mathematical physics. 2015. Vol. 55, no. 7. P. 1077-1093.

19. Калиткин Н. Н. Численные методы. СПб. : БХВ-Петербург, 2011. 592 с.

20. Korytov M.S., Breus I.V. Development of the mathematical description of a overhead crane in large spatial movements, taking into account the dissipation of swing energy // International Journal of Mechanics and Control. 2018. Vol. 19. No. 02. Pp. 39-44.

УДК 621.873.1

МОДЕЛЬ СФЕРИЧЕСКОГО МАЯТНИКА С ПОДВИЖНОЙ ТОЧКОЙ ПОДВЕСА В ЗАДАЧЕ ПРОСТРАНСТВЕННОГО ПЕРЕМЕЩЕНИЯ ГРУЗА ГРУЗОПОДЪЕМНЫМ КРАНОМ

ПРИ ОГРАНИЧЕНИИ КОЛЕБАНИЙ

MODEL OF THE SPHERICAL PENDULUM WITH A MOBILE POINT OF A SUSPENSION IN THE TASK OF THE SPATIAL MOVEMENT OF A CARGO OF A LOADING CRANE

WITH LIMITATIONS OF SWAYS

М. С. Корытов1, В. С. Щербаков1, В. В. Титенко2, В. Е. Беляков3

'Сибирский государственный автомобильно-дорожный университет, г. Омск, Россия 2Омский государственный технический университет, г. Омск, Россия 3Военная академия материально-технического обеспечения имени генерала армии А. В. Хрулева Министерства

обороны Российской Федерации, г. Омск, Россия

M. S. Korytov1, V. S. Shcherbakov1, V. V. Titenko2, V. E. Belyakov3

' Siberian automobile and highway university, Omsk, Russia 2 Omsk State Technical University, Omsk, Russia 3 Military Educational Institution of Logistics named after General of the Army A. V.Khrulyov of the Ministry of Defense

of the Russian Federation, Omsk, Russia

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

Ключевые слова: грузоподъемный кран, груз, ограничение колебаний, раскачивание, сферический маятник.

DOI: 10.25206/2310-9793-7-1-104-110

I. Введение

Грузоподъемные краны (ГК) широко используются в качестве механизмов перемещения грузов во многих строительных и промышленных областях [1]. Чтобы повысить эффективность и безопасность ГК, возникающее при их движении, раскачивание груза на гибком канатном подвесе должно быть ограничено [2].

Проблема ограничения раскачивания грузов, перемещаемых ГК, была широко исследована. В работе [3] были предложены аналитические зависимости оптимального и квазиоптимального режимов управления маятниковой системой, имеющей одну угловую степень свободы с подвижной точкой подвеса для задачи наискорейшего разгона (торможения) с гашением колебаний. На скорость и ускорение точки подвеса наложены ограничения. Недостатки известного способа - это рассмотрение малых колебаний маятника вокруг положения равновесия, неопределенность предельного значения угла отклонения грузового каната ГК от вертикали. Оптимальное управление носит релейный характер: ускорение точки подвеса принимает граничные и нулевое значения. Гашение колебаний маятника происходит не на всем интервале времени рабочего цикла, а лишь к концу разгона или перемещения системы [3, 4].

Применение для решения указанной проблемы ПД- и ПИД-регуляторов [5, 6], аппарата нечеткой логики и нейронных сетей [7, 8] и shaping-алгоритмов [9, 10] при всех различиях в подходах создает сравнительно большую погрешность реализации как угла отклонения грузового каната ГК от вертикали, так и линейных координат перемещения груза. Неуправляемая компонента маятниковых колебаний груза подавляется не полностью. Рассматриваются в основном системы с одной угловой степенью свободы [11]. Время перемещения при гашении колебаний, как правило, увеличивается.

Целесообразна разработка алгоритма пространственного перемещения груза по криволинейной траектории, который при заданных ограничениях в виде максимальных скорости и ускорения подвижной точки подвеса груза на ГК, который синтезировал бы непрерывное (бесступенчатое, не релейное) управление точкой подвеса при помощи частотно-регулируемых приводов ГК. Такими приводами оснащается ряд изготовляемых в настоящее время ГК. Алгоритм должен синтезировать траекторию с заданными углами начальных и конечных перемещений груза в неподвижной системе координат в плане, а также учитывать возможность больших углов отклонения грузового каната от гравитационной вертикали, что позволит максимально повысить скорость перемещения и производительность ГК. Целесообразно использовать в качестве математического описания груза, перемещаемого ГК на гибком канатном подвесе, дифференциальные уравнения в сферических координатах. Они подходят для описания больших угловых перемещений грузов различными видами ГК [12]

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

Для возможности моделирования пространственных перемещений груза ГК, описываемых большими углами отклонений грузового каната от гравитационной вертикали (свыше 5°), использовалась система известных нелинейных дифференциальных уравнений второго порядка вида [13, 14]

róe = (Z/ 2)гоф sin(29)-(l/Z)(g - sin8 + (3q coscp + 3r2 sincp)-cos8);

0 = Год;

(1)

соф = -2гоегофс1|>9 +((3г! sincp -х2 coscp)/(Z-sin8)); ф = гоф,

где L - длина подвеса грузового каната от точки подвеса до центра массы груза; 0 - угол отклонения груза и каната от гравитационной вертикали; ф - угол поворота маятниковой системы в горизонтальной плоскости; 8 = сое,ф = озф - обобщенные скорости; g - ускорение свободного падения; х,. х2 - ускорения точки подвеса локальной вдоль двух горизонтальных осей X1 и X2 неподвижной прямоугольной системы координат соответственно (рис. 1).

При исследованиях были приняты допущения об отсутствии ограничений на амплитуду угла 0, о постоянстве длины грузового каната L в процессе перемещения груза, о бесступенчатом характере регулирования скоростей хъ х2 и ускорений х,. х2 разгона и торможения точки подвеса груза (использование частотно-регулируемого электропривода) вдоль осей X1 и X2 и о пренебрежимо малом влиянии массы перемещаемого груза и подвижных звеньев ГК на управляемые параметры скоростей и ускорений точки подвеса.

В качестве примера, имеющего широкое практическое применение, рассматривалась траектория, начальное направление которой параллельно движению ГК вдоль оси X1, а конечное параллельно движению ГК вдоль оси X2.

В качестве исходных данных задачи помимо перечисленных выше постоянных параметров L и g выступали конечные требуемые горизонтальные координаты груза x1end и x2end; начальный угол поворота маятника в горизонтальной плоскости q>start: конечный угол поворота маятника в горизонтальной плоскости ф„,и/: максимально допустимые скорости х, ]п11 , х2нт и ускорения xllim , x21im точки подвеса груза.

Груз из состояния покоя на вертикальном канатном подвесе должен быть перемещен ГК на заданные расстояния x1end и x2end вдоль осей X1 и X2 за время T. После перемещения (в момент времени T) груз также находится в точке с координатами [x1end; x2end] в состоянии, близком к состоянию покоя (вертикальное положение грузового каната, отсутствие остаточных колебаний).

Приняты обозначения вспомогательных постоянных параметров алгоритма: абсолютных пороговых (ниже которых краевые условия считаются выполненными) значений угла каната относительно вертикали, его скорости и ускорения (),. 0, и 0, соответственно, абсолютных пороговых значений отклонения координаты точки подвеса вдоль каждой из осей A'i иЛ'2, ее скорости и ускорения л.-,. xt и х! соответственно.

Целевая функция (2), краевые условия (3) и ограничения (4) задачи имеют вид, приведенный ниже.

T ^ min. (2)

Вид целевой функции (2) обусловлен необходимостью повышения производительности ГК за счет минимизации времени перемещения груза.

lei п <е,. |ё| <ё,, |ё| <ё,,

I 1/=0, Г ' | |/=0, Г I \t=0,T '

lXl|/=0 -Х/' lxl ~Xlend\t=T -Х/' lXl|/=0. Г -Xi' lXl|/=0. Г

lX2|i=o -Xf |Х2 ~X2end\t=T ~Xt> |Х2^=0 - Xf |Х21/=0. Г ~Xt •

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

1*1 (0| ^ *111т, 1*1 (0| ^ *111т, 1*2 (0| ^ *2Нт, |*2 (0| ^ *2Нт V Г £ [0, Г]- (4)

Ограничения следуют из невозможности превышения максимально допустимых скоростей и ускорений подвижных звеньев ГК.

III. Теория

Временные зависимости ускорений точки подвеса груза х = f{t) были получены в аналитическом виде дифференцированием предложенных временных зависимостей 0=ДО угла отклонения груза и каната от гравитационной вертикали и ф=Д0 угла поворота маятниковой системы в горизонтальной плоскости [15, 16]:

9(0 = A/(e^1 +1) + AI(e"^2 +1);

(5)

Ф(0 = Фе„, /(e^-j +1);

(6)

где А - амплитуда угла 0 отклонения груза и каната от гравитационной вертикали во время перемещения груза, рад; k0 - коэффициент крутизны нарастания и спада значения угла 0; кф - коэффициент крутизны нарастания значения угла ф; cj0 и c20 - временные значения локальных центров нарастания и спада угла 0 грузового каната ГК, с; сф - временное значение центра нарастания угла ф грузового каната ГК, с.

Временные зависимости углов 0=ДО и Ф=Д0 представлены в виде суммы элементарных сигмоидальных функций. Причем выражение (5) стремится к нулю на минус и на плюс бесконечности аргумента времени, а выражение (6) - только на минус бесконечности. На плюс бесконечности функция (6) стремится к ф„,и/. Стремление (5) к нулю обеспечивает при задании достаточно малых значений О,. 0, и 0, отсутствие существенных колебаний груза как в начальный, так и в конечный моменты рассматриваемого перемещения.

Выражения (5) и (6) могут быть дифференцированы в аналитическом виде, что позволило получить выражения скоростей и ускорений углов наклона каната 0 и ф:

0(0 =

(7)

0(0 =

Акп2 екв(с-в~п Ак.2 екв(С1в~п 2 Акп2 е2кв(С1в~° 2 Ак<? е2кв(с-в~п

- + -

(еАе(с2е-0 (eA0(ci0_í)+1)2 (eA"e(Cl0_í)+1)3 (екв(с2в-*)

(8)

1)2

(9)

ф(0 =

2(Vendk<?~e

2„2К (Сф-t)

ФеМкф e

2 „кФ (сф-)

(екф(сф f) +i)3

(екф(сф f) +i)2

(10)

Выражения скоростей и ускорений углов (7)—(10). в свою очередь, позволяют вывести из (1) аналитические зависимости ускорений точки подвеса х,. х-, от остальных переменных системы уравнений (1). Они имеют вид:

( L ф2 sin(2 6) sin(0) сов(ф) - 2 Z 0 sin(0) сов(ф) - ^ -2 g (sin(8))2 со^ф) + 2 L ф cos(8) (sin(0))2 вЦф) + +2 L 0 ф sin(2 0) cos(0) вЦф)

2 cos(0) sin(0) (( cos^))2 + (sin^))2 )

(11)

2 g (sin(0))2 8ш(ф) + 2 L 0 sin(0) 8ш(ф) -—L ф2 sin(20)sin(0)sin^) + 2.Z^cos(0)(sin(0))2 соз(ф) + +2 L 0 ф sin(2 0) cos(0) соз(ф)

2 cos(0) sin(0) |(соз(ф))2 + (зЦф))2 j

(12)

2

Подстановкой выражений (5)—(10) в (11) и (12), ускорения подвеса х,. х-, могут быть выражены в развернутом виде исключительно через константы и переменную времени t. Ввиду громоздкости указанные развернутые выражения не приводятся.

Наконец, известные временные зависимости ускорений (11), (12), скоростей и перемещений точки подвеса позволят перемещать подвижные звенья ГУ (где находится точка подвеса) по траекториям, обеспечивающим отсутствие колебаний груза в начальный и в конечный моменты времени.

Вывод аналитических выражений интегралов х, (I). х1(г) и х2(/). х2(г) затруднен. Поэтому векторы дискретных значений ускорений точки подвеса хъ х,, полученные по зависимостям (11), (12) путем подстановки численных значений, полученных по формулам (5)-(10) для различных моментов времени t е[0, Т] с определенным шагом дискретизации может быть дважды проинтегрирован при помощи известного численного метода трапеций [16].

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

ч (1 /Л (7) Л/ ч (1 / сН (8) к, ч _ •• , ч ¡Л . , ч \Л

/ч d/dt (9) ,/ч d/dt (10) ,./ч .. ^ ¡dt . ^ ¡dt (p(t)--— __>. x2(?) J > х2 it- ) —-->x2(t) (13)

Аналитические выражения Численное интегрирование

Начальная часть схемы (13) вычисляется по аналитическим зависимостям, поэтому получаемые значения перемещений точки подвеса Xj(t), x2 (t) характеризуются высокой точностью. Схема (13) также характеризуется высоким быстродействием и вычислительной устойчивостью.

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

Д/0 - интервал времени от центра элементарной возрастающей (убывающей) сигмоиды изменения угла 9 до момента времени ранее (позже), когда значение функции составляет достаточно малую долю (1-P) от амплитуды угла;

Дф - интервал времени от центра элементарной возрастающей сигмоиды изменения угла ф до момента времени ранее, когда значение функции составляет достаточно малую долю (1-P) от амплитуды угла;

Д/02=С20-С10 - интервал времени между центрами двух элементарных сигмоид изменения угла 0 на рост и на спад.

Коэффициенты крутизны нарастания и спада значений углов 0 и ф в элементарных сигмоидах (5), (6) связаны со значением P зависимостями

к, = 1п(1/р - 1)/(-д/0); к = 1п(1/р - 1)/(-дд.

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

В качестве дополнительного условия, уменьшающего размерность решаемой задачи, было принято условие совпадение временных значений центра сф единственной элементарной сигмоиды угла ф и общего центра (с10+с20)/2 двух элементарных сигмоид угла 0 на рост и на спад:

сф=(с10+с20)/2.

IV. Результаты экспериментов

Было установлено, что изменение положительных значений каждого из трех временных параметров Д/0, Д/"ф и Ai02 оказывает неоднозначное влияние на максимальные ускорения, максимальные скорости и конечные перемещения точки подвеса груза xllim, x21im - *iiim - *2iim - xiend - x2end ПРИ различных сочетаниях значений Дi0, Дф и Д/02. Поэтому для выполнения всех краевых условий (3) и ограничений (4) задачи был применен Симплекс-метод оптимизации функции четырех аргументов [A; Д0; Дф; Д/02] в сочетании с методом штрафных функций [17].

Значение минимизируемой целевой функции определялось выражением

S = |Х1 ~ Х1 end I + |Х2 (T) ~ X2end | + |*2 ~ X21im | + ^sh ' Т .

где ksh - эмпирический коэффициент.

Штрафные функции вычислялись согласно выражениям:

Л = 1*1 (Т) - *Шт I • еСЛИ 1*1 (Г)| > *Шт • Л = 1*1 (Т) " *11пп |" Кк еСЛИ |*1 (Г)| > *]

Л = |*2 *21пп | * еСЛИ ^Ъ^)^*

1Нт •

21пп '

Полное выражение целевой функции с учетом штрафных функций имело вид:

/(А; ; ; ^е2) = £ + / + /2 + / .

(14)

Перечисленные выше постоянные исходные параметры задачи в рассматриваемом примере принимали значения Ь =10 м; g =9,81 м/с2; х1икр10 м; х2е„а=8 м; ф^=0°; фик/=180°; х1тах=1,5 м/с; х2Ьп=1,2 м/с;

1Нт

=0,4 м/с2, Зг2Цт =0,3 м/с2; к,ь =0,1, Р=0,999.

1 0 -1 -2 -3

10

8. град; 9. град/с; 9. град/с

200

100

0 5 10 15 /, с 25

а)

„2

-100

ф, град; ф , град/с; ф, град/с2

хь м; £ м/с; х м/с

\ Х1

Х1 /

А, ^ - - — -

х2, м

0 5 10 15 /, с 25

б)

х2, м; х2, м/с; х2, м/с

Х1 у/

1

г

5 10 15 с 25 0 5 10 15 с 25

б) г)

4 д) 6

8 х:, м 10

Ф

0

5

0

0

8

0

2

Рис. 2. Пример временных зависимостей функций изменения углов 9(/) и ф(/) грузового каната, их первых двух производных и соответствующих им линейных перемещений точки подвеса и груза в пространстве

Начальные значения аргументов при оптимизации принимали значения [4=1°; Д/е=5 с; Дф=5 с; Д/ег=5 с]. На рис. 2 приведены графические результаты моделирования оптимального по (14) перемещения ГК и груза: а - временные зависимости угловой координаты е и ее первых двух производных; б - временные зависимости угловой координаты ф и ее первых двух производных; в - временные зависимости линейного перемещения моста и его первых двух производных; г - временные зависимости линейного перемещения грузовой

тележки и его первых двух производных; д - вид в плане перемещений точки подвеса груза на рабочей площадке ГК (за счет перемещений ГК) и самого груза.

VI. Выводы и заключение

Полное время перемещения в рассматриваемом примере составило 7=24.92 с. Максимальные достигнутые во время перемещения ускорения и скорости составили jrlmax =0,819 м/с2, х2тлх =0.973 м/с2, х1тах =1,416 м/с, х2тах =0,642 м/с. Все они меньше соответствующих предельных значений jfllim, x21im, xllim, x21im.

Найденные оптимальные значения временных параметров равны Д/е=8,046 с, Д/ф=2,983 с, Д/е2=8,835 с. Оптимальное значение амплитуды отклонения угла от гравитационной вертикали A=2,36 град.

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

Список литературы

1. Inomata A., Noda Y. Fast trajectory planning by design of initial trajectory in overhead traveling crane with considering obstacle avoidance and load vibration suppression // Journal of Physics: Conference Series. 2016. Vol. 744. No. 1. P. 012070.

2. Matsunaga M., Nakamoto M., Yamamoto T. A sound-based measurement of sway angle for anti-sway control of overhead crane // Journal of robotics networking and artificial life. 2018, Vol. 4, no. 4. P. 322-325.

3. Черноусько Ф. Л., Акуленко Л. Д., Соколов Б. Н. Управление колебаниями. М.: Наука, 1980. 383 с.

4. Reshmin S. A., Chernous'ko F. L. A time-optimal control synthesis for a nonlinear pendulum // Journal of Computer and Systems Sciences International. 2007. Vol. 46, no. 1. P. 9-18.

5. Mohd Tumari M. Z., Shabudin L., Zawawi M. A., Ahmad Shah L. H. Active sway control of a gantry crane using hybrid input shaping and PID control schemes // IOP Conference Series: Materials Science and Engineering. 2013. Vol. 50(1). P. 012029. DOI: 10.1088/1757-899x/50/1/012029.

6. Korytov M. S., Shcherbakov V. S., Titenko V.V . Comparative analysis of methods of cargo vibration damping moved by overhead crane // Journal of Physics: Conference Series. Vol. 1050, no 1. 2018. P. 012038. DOI: 10.1088/1742-6596/1050/1/012038.

7. Shehu M. A., Li A., Huang В., Wang Y., Liu B. Comparative Analysis of Neural -Network and Fuzzy AutoTuning Sliding Mode Controls for Overhead Cranes under Payload and Cable Variations // Journal of control science and engineering. 2019. P. 1480732. DOI: https://doi.org/10.1155/2019/1480732.

8. Huang X., Gao H., Ralescu A.L., Huang H. Adaptive hierarchical sliding mode control based on fuzzy neural network for an underactuated system // Advances in mechanical engineering. 2018. Vol. 10, no 9. P. 1687814018799554. DOI: 10.1177/1687814018799554

9. La V. D., Nguyen K. T. Combination of input shaping and radial spring-damper to reduce tridirectional vibration of crane payload // Mechanical systems and signal processing. 2019. Vol. 116. P. 310-321. DOI: 10.1016/j.ymssp.2018.06.056.

10. Blackburn D., [et al.]. Command Shaping for Nonlinear Crane Dynamics // Journal of Vibration and Control. 2010. N 16. P. 477-501.

11. Korytov M. S., Shcherbakov V. S., Titenko V. V. Analytical solution of the problem of acceleration of cargo by a bridge crane with constant acceleration at elimination of swings of a cargo rope // Journal of Physics: Conference Series. 2018. Vol. 944. no. 1. P. 012062. DOI: 10.1088/1742-6596/944/1/012062

12. Samin R. E., Mohamed Z. Comparative Assessment of Anti-Sway Control Strategy for Tower Crane System // AIP Conference Proceedings. 2017. Vol. 1883. UNSP 020035. DOI: 10.1063/1.5002053.

13. Блехман И. И. Вибрационная механика. М.: Физматлит, 1994. 400 с.

14. Бутиков Е. И. Необычное поведение маятника при синусоидальном внешнем воздействии // Компьютерные инструменты в образовании. 2008. № 2. С. 24-36.

15. Mitchell T. M. Machine Learning. WCB/McGraw-Hill, 1997. 414 p.

16. Крылов В. И. Приближенное вычисление интегралов. М.: Наука, 1967. 500 с.

17. Васильев Ф. П. Методы оптимизации. М. : Факториал Пресс, 2002. 824 с.