Исследование методов параметрической интерполяции для сглаживания кусочно-линейной траектории движения инструмента Текст научной статьи по специальности «Физика»

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ И

ПРОЦЕССЫ УПРАВЛЕНИЯ N. 4, 2021 Электронный журнал, рег. Эл № ФС77-39410 от 15.04.2010 ISSN 1817-2172

http://diffjournal. spbu. ru/ e-mail: _ jodiff@mail.ru

Общая теория управления

Исследование методов параметрической интерполяции для сглаживания кусочно-линейной траектории движения

инструмента

1 SH 1 Jjiiji 1 *** 1

Зеленский А.А. , , Абдуллин Т.Х. , , Дубовсков В.В. , , Илюхин Ю.В. , v

1 ФГБОУ ВО "МГТУ "СТАНКИН"

*zelenskyaa@gmail.com **everestultimate@yandex.ru ***dtnt1121@gmail.com *** *ilyv-2@mail.ru

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

Ключевые слова: параметрическая интерполяция, алгоритм Рунге-Кутты, разложение в ряд Тейлора, сглаживание траектории, контроллер движения, высокоскоростная обработка.

1. Введение

Параметрические кривые нашли широкое применение в системах автоматизированного проектирования (САПР) при проектировании деталей, описываемых поверхностями свободной формы [1]. Основная проблема обработки таких изделий заключается в том, что траектория перемещения режущего инструмента в CAM-системе (Computer-aided manufacturing) представлена в параметрической форме и регулируется стандартами STEP (STandardfor Exchange of Product model data) или IGES (Initial Graphics Exchange Specification) [2]. В свою очередь, большинство систем управления не соответствуют этим стандартам, а CAM-система в таком случае генерирует либо линейные (G01), либо круговые кадры (G02, G03), описываемые стандартом ISO6983 [3]. В результате этого инструмент будет перемещаться вдоль кусочно-линейной траектории с геометрической непрерывностью G0 (два сопрягаемых линейных участка имеют одинаковые координаты в граничной точке).

Отработка таких участков траектории по алгоритму линейной интерполяции в случае мгновенного изменения направления движения с заданной скоростью подачи будет сопровождаться резкими скачками контурного ускорения. В результате будут происходить динамические удары по узлам станка или промышленного робота [4,5], ухудшится качество механообработки и возможно проявление ограничений, действующих в позиционном контуре сервоусилителей. А это может привести к аварийному останову одного из звеньев механизма. Поэтому в системе управления необходимо учитывать кинематические ограничения каждого из звеньев механизма на допустимый рывок и ускорение. Однако тогда отработка кусочно-линейной траектории будет выполняться лишь с незначительными скоростями подач, и это существенно снизит производительность технологического оборудования [4-7].

Стремление к повышению точности аппроксимации поверхности изделия в CAM-системе путём увеличения количества кадров в управляющей программе только усугубляет проблему, т.к. длина каждого отрабатываемого сегмента траектории уменьшается. Это приводит к увеличению размера передаваемых данных, повышению нагрузки на аппаратную часть системы управления и потребует большего объёма оперативного и постоянного запоминающего устройств [6].

Решение обозначенных проблем заключается в использовании алгоритма сглаживания траектории инструмента с заданным допуском построения сплайна, при котором исходная траектория будет преобразовываться в линейно-сплайновые участки. Такой метод требует реализации параметрической интерполяции в системе управления, где кинематические ограничения будут связаны с кривизной сплайна [8-11]. Это способствует повышению контурной скорости подачи на сопрягаемых участках и производительности обработки [12,13].

Параметрическая интерполяция сопровождается расчётом соответствующего параметра сплайна u в каждый дискретный промежуток времени, а кривая представляет собой функцию С (и) = [x(u),y(u),z(u)], где ие [a, b]. В качестве сопрягающей кривой обычно используют B-сплайны [10,11,14], кривые Безье третьего [15] и пятого порядков [16], кривые годографа Пифагора [17], либо неоднородные рациональные B-сплайны (кривая NURBSa) [18]. Среди этих типов параметрических кривых отличительной особенностью обладают кривые годографа Пифагора, в которых длина дуги может рассчитываться аналитически исходя из заданных контрольных точек сплайна и некоторого полинома &(и), связанного с первыми производными годографа кривой С'(и) = [x(u),y(u),z(u)] с помощью уравнения

х(и)2 + у(и)2 + ¿(и)2 = а(и)2, (1)

где а (и) = |C'(u)| = ds/du - определяет параметрическую скорость кривой, ds - абсолютное приращение пути вдоль дуги кривой, du - соответствующее изменение параметра. Тем не менее,

многие CAM/CAD системы не поддерживает данные кривые [4], а в работах, где рассматриваются открытые, либо экспериментальные системы управления, по большей части применяют кривые Безье, B-сплайны или кривые NURBSа. В закрытых проприетарных системах (SIEMENS 840D, FANUC) также реализуется функция локального сглаживания, однако применяемый тип кривых не раскрывается [19,20].

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

Одним из вариантов решения проблемы является применение аппаратных ускорителей на базе ПЛИС и построение на их основе распределённых вычислительных систем [24,25].

Другая проблема состоит в том, что дискретное приращение пути вдоль кривой и соответствующее приращение параметра сплайна характеризуется нелинейной зависимостью [26]. При этом явная аналитическая связь между параметром сплайна и его длиной отсутствует. Это приводит к значительным погрешностям при расчёте новой интерполируемой точки, а значения скоростей подач подвержены колебаниям, которые негативным образом сказываются на качестве обрабатываемой поверхности. Применение аппроксимирующих и итерационных методов поиска параметра сплайна [27-31], либо параметризация длины дуги [26,32] полиномами высших порядков частично решают данную проблему, однако не происходит полного устранения флуктуации значений скоростей подач. При параметризации длины дуги полином должен иметь по крайней мере пятую степень, которая позволяет наложить граничные условия нулевого, первого (ds/du) и второго порядков (d2u/ds2) в начальной и конечной точках сплайновых сегментов, обеспечивая тем самым непрерывность их стыковки. Дальнейшее повышение степени полинома даёт возможность точнее выполнить его приближение с исходной геометрией заданной кривой и уменьшить флуктуации скоростей подач. Таким образом, задача сглаживания траекторий остаётся актуальной.

Стоит также отметить, что в алгоритме параметрической интерполяции должна контролироваться хордовая ошибка построения сплайна 8 (рис. 1а), т.к. большая скорость подачи будет приводить к меньшей точности обработки контура [6], а наибольшее значение хордовой ошибки будет наблюдаться в сегментах, обладающих меньшим радиусом кривизны.

Применяемая схема сглаживания кусочно-линейной траектории непосредственно влияет на качество и эффективность обработки контура. Например, вставленный кубический B-сплайн с тремя контрольными точками обеспечивает пропорциональность первых производных по параметру двух сопрягаемых сегментов и равное направление касательных векторов в общей точке сопряжения линейного и параметрического участков, т.е. выполняется условие геометрической непрерывности G1 [33]. Кроме того, в данном случае не контролируется ограничение на нормальную составляющую рывка [10], а в линейно-сплайновом переходе будет возникать скачкообразное изменение контурного ускорения и скорости. Для устранения данного недостатка сопряжение линейно-сплайновых сегментов должно выполняться как минимум с геометрической непрерывностью G , где первая и вторая параметрические производные на границе двух сопрягаемых участков пропорциональны между собой. В этом случае соблюдается равенство кривизны и направления касательных векторов в общей точке соединения обоих сегментов. При этом сама переходная параметрическая кривая должна иметь непрерывность C2 вдоль всего сглаживаемого участка, а выбранный алгоритм параметрической интерполяции должен обеспечивать высокую точность расчёта параметра сплайна.

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

кривая представлена в виде симметричного кубического В-сплайна с пятью контрольными точками [10,11]. Проведен краткий обзор и сравнительный анализ наиболее распространённых методов параметрической интерполяции, обладающих разной точностью аппроксимации параметра кривой.

Все вычислительные эксперименты учитывают алгоритм планирования подач [34], за формирование которого отвечает ряд следующих взаимосвязанных программных компонентов: геометрический модуль сглаживания траектории, двунаправленный алгоритм предпросмотра кадров и отдельный модуль s-образного планирования подач, включающий оптимизацию скорости. Описание перечисленных алгоритмов выходит за рамки данной статьи.

2 Обзор наиболее распространённых методов параметрической интерполяции

Основными входными данными для интерполятора, содержащего алгоритмы линейной и параметрической интерполяции, являются дискретное приращение пути Лб^, определенное на этапе Б-образного планирования подач, координаты контрольных точек, суммарный пройденный путь и длины линейного и параметрических участков.

Обобщённая формула параметрической интерполяции вдоль кривой, соответствующая пройденному пути ЛЬу, имеет следующий вид:

= иу + Ли.], (2)

где значение параметра последующей интерполируемой точки С(и^+г), и^ - текущее

значение, соответствующая точке С(и^), Ли- приращение параметра за дискретный период Т3. Для оценки хордовой ошибки построения сплайна используются некоторые упрощения с построением вспомогательной окружности с радиусом . Тогда последующая интерполируемая точка Р(и]+1) будет располагаться на дуге (рис. 1б) и, учитывая пройденный путь ЛЬ], возможно рассчитать приближенное значение 8у, соответствующая отрезку АВ [8] по формуле

= Ъ - №

~2 (3)

(ЛЬ^ 2 v 7

Кривизна сплайна К] и радиус кривизны в точке С (и.]) рассчитываются по формуле

д ГМГ , 1 (4)

"К'ММГМИ ' я/

где С'(щ), С '(и) первая и вторая производные функции параметрической кривой С (и), ||. || -евклидова норма. Контроль хордовой ошибки построения сплайна обеспечивается на предыдущем этапе с использованием уравнения (3) при помощи двунаправленного алгоритма предпросмотра кадров. При этом выполняется оценка допустимой скорости подачи в критической точке сплайна, соответствующей максимальной кривизне ктах:

2 I 28 (5)

^minк = -п I __ З2,

Тя л Ктах

где 8 - заданная в системе управления хордовая ошибка построения сплайна.

Недостаточно точный расчёт параметра является причиной возникновения колебаний скоростей подач, при которых реальное смещение ЛЬу вдоль кривой отличается от запланированного дискретного смещения Лб^ . Флуктуации скорости наиболее заметны в самом простом случае при естественной интерполяции, когда параметр и имеет пропорциональную зависимость от желаемого смещения [26]:

Н

и = —s,

L

(6)

где Ь - длина всего сплайна, 5 - смещение вдоль дуги кривой, Н - максимальное значение параметра сплайна. Данный тип интерполяции применяется только для кривых, описываемых полиномиальными уравнениями, при сопряжении которых не будут выполняться граничные условия по первой производной.

Рис. 1 Обобщенная схема параметрической интерполяции (а), схема расчёта хордовой

ошибки Sj (б).

При расчёте параметра сплайна удовлетворительные результаты получаются только в случае, если зависимость длины дуги от параметра сплайна близка к линейной. Однако было доказано [35], что математически невозможно сгенерировать сплайн с полностью линейной зависимостью, иначе он будет представлять собой прямую линию с отсутствием кривизны. Относительные колебания скоростей подач для данного метода могут достигать десятка процентов [32]. Но для высокоскоростной обработки их допустимые значения составляют 0.001...0.1% [29-31]. Как показывает опыт предыдущих работ, добиться такого результата возможно за счёт применения аппроксимирующих методов, обладающих разной точностью приближения параметра сплайна.

Среди методов, показывающих приемлемые результаты, стоит выделить параметризацию длины дуги сплайна корректирующим полиномом подачи (англ. feed correction polynomial) [26,32], обеспечивающим корректное сопоставление желаемого смещения с параметром сплайна. Интерполяция реализуется вдоль кривой NURBS^ а его длина на интервале [а, Ь] при а <и <Ъ будет выражаться через определённый интеграл:

l(u) = fb ^х(и)2 + у(и)2 + z(u)2,

(7)

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

1(а,Ъ) = (h/3)(f(d) + 4f(c)+f(b)),

(8)

где с = (а + Ь)/2 - середина интервала [а,Ь], к = (Ь — а) - длина интервала, /(и) - первая производная функции длины, рассчитываемая по формуле

/(и) = 1'(и) = ^х(и)2 + у(и)2 + г(и)2. (9)

Используя метод деления отрезка пополам и ограничиваясь заданным допуском аппроксимации е, получим два новых подинтервала [а1,Ь1] и [а2,Ь2], длина которых также рассчитывается по соотношению (8). Прекращение итерационного процесса разбиения произойдет в том случае, если будет выполняться следующее условие Ща^Ь^) + 1(а2,Ь2) — 1(а,Ь)^10 < £. В противном случае, разбиение и расчёт длины интервала продолжится до удовлетворения заданного допуска. В результате формируется массив из М строк, который включает в себя длины каждого из подинтервалов 1(иу), параметр сплайна и суммарное смещение (41,51). Полученный массив позволяет сгенерировать функцию и = /(з), которая служит эталоном по подбору корректирующего полинома подачи. Искомый многочлен может иметь как седьмой, так и девятый порядок и описывается следующим соотношением (рис. 2):

и = /(Б) = А959 + Л858+... +А0,

(10)

где 0 < 5 < Ь, А9 ...А0 - подобранные коэффициенты, и = [й0,...,йм]т- рассчитываемые параметры сплайна. После расчёта полинома новые значения его точек (щ, ${) сравниваются с исходными данными через среднеквадратичную ошибку. В случае, если она не будет удовлетворять заданному допуску £а, то весь набор точек делится пополам, и для двух новых интервалов расчёт многочлена запускается повторно. Среднеквадратичная ошибка рассчитывается по формуле [22]:

N.

а

I

1 = 0

(щ — щУ

К

(11)

< е„

где Ы3 < М - количество точек после разбиения интервала.

х Действительное смещение I-

Смещение при помощи корректирующего полинома подачи Рис. 2 Применение корректирующего полинома подачи.

Приведенные в работе [32] эксперименты показывают, что максимальное относительное отклонение скорости подачи от запланированного значения для тестового веерообразного контура составляет порядка 0.13%. Между тем данный подход требует довольно затратной предварительной подготовки с использованием методов численной оптимизации. В частности, каждая итерация по нахождению многочлена представляет собой оптимизационную задачу квадратичного программирования с линейными ограничениями, решаемую при помощи метода множителей Лагранжа. Вследствие этого на практике данный метод реализован в постпроцессорах

CAM/CAD систем. Корректирующий полином встраивается в кадры управляющей программы, по которому затем осуществляется аппроксимация параметра сплайна в реальном времени.

Таким образом, немаловажным свойством применяемых алгоритмов является исключение предварительной подготовки и возможность реализации алгоритма в режиме реального времени непосредственно в самой системе управления. Среди таких способов стоит выделить разложение в ряд Тейлора второго порядка [10,27] и классический метод Рунге-Кутты (РК) четвертого порядка [28,29], относящиеся к прямым способам приближения или методам преобразования без обратных связей.

Из-за усечения членов высших порядков в ряде Тейлора при отработке траектории с большой кривизной будут возникать нежелательные колебания в скоростях подач (до 0.135%) [30]. Соответствующее выражение, которое связывает последующий параметр сплайна Uj+1 с его текущим значением Uj, записывается следующим образом:

+ т • + (Ts2Y + (12)

Uj+1 = Uj + TsUj + JUj + ч.в. п., где Uj, Uj раскрываются через формулы:

dU

U = — d

V/ V/ V/ (13)

= IIГ

t=tj \\C'(Uj)W Jx2 + y2 + z2

du

где x = dCx(U)/dUlu=u ,y = dCy(ü)/dU\ , y = dCy(U)/dU\ ;

u= u u= u

d2 U

U = ~—r d 2

= ]

( V^) -d2s/dU2 _ (14)

(ds/dti)3

2

= (У/) (х-х + у-у + 1-2) = 2 (С'(ц), С"(и])) ^х2+у2 + г2 (]) \\СХи)\\4 '

где х = d2Сx(и)/dи2\u=U|., у = d2Сy(и)/dи2\ г = d2Сz(и)/dи2\u=u. .

' И_ uj J

С учетом соотношений (13) и (14) уравнение (12) записывают с использованием дискретных приращений [10,28]:

1 (х-х+у-у + г-ё)2 (15)

= и + ^2+у2+г2Л51 - 2-(х2+у2+г2)2 ^ '

где АЗ] = У^Т3, V'¡* = (У}+1 + У])/2 - действующая скорость подачи на дискретном шаге Т5, У]+1 и У - запланированные значения скоростей подач в интерполируемых точках. Формулы для расчёта колебаний скоростей подач приобретают вид

^ = ^М х 100%, V] = П^-чиЯ = *Л (16)

* Б * Б

где V] — реальная скорость подачи на выходе интерполятора.

Как показано в работе [30], аппроксимация параметра рядами Тейлора была улучшена при помощи схемы компенсации длины дуги, а время расчёта при этом увеличилось незначительно.

Отклонения скоростей подач укладываются в приемлемые значения до 0.065%. В основе метода лежит положение о том, что флуктуации скоростей подач вызываются расхождением реальной и идеальной траекторий сплайна. Из соотношений (16) можно заключить, что полная компенсация колебаний возможна при выполнении следующего равенства:

\\С(и]+1) - С(и])\\ = ЛЬ= Щ = Л8}. (17)

Однако в действительности оно выполняться не будет, т.к. при параметрической интерполяции пройденный путь отрабатывается вдоль отрезка, а не точно по контуру дуги сплайна. Как можно заметить, реальная скорость подачи V] вдоль отрезка имеет меньшее значение, чем скорость вдоль дуги. Это показывают формулы

АВ

Т

>

\\АВ\\

(18)

Т

V > у

где символом У У обозначена длина дуги. Поэтому для компенсации длины дуги сплайна последующую точку интерполяции сдвигают вдоль кривой (рис. 3а), т.е. \\АС\\ = У]*Т5. Далее длина дуги сплайна локально аппроксимируется окружностью кривизны от текущей точки интерполяции. Полагая, что длины отрезков АС и АБ равны между собой, длина дуги окружности

АБ может использоваться для аппроксимации длины дуги сплайна АС (рис. 3б). Исходя из этого, новое соотношение для дискретного смещения выглядит следующим образом:

Л§] =

А С

АБ

У?Т3

= Я]в = 2Я]агсв\п

2 К]

у;т +

(у/Т)3

2\Щ

(19)

г\

где радиус кривизны рассчитывается из соотношения (4). Оно подставляется в выражение (15).

Рис. 3 Схема компенсации длины дуги: (а) коррекция заданного положения, (б) приближенная оценка длины дуги с использованием окружности кривизны; О'- центр окружности, К и в - радиус кривизны и центральный угол соответственно.

По сравнению с разложением в ряды Тейлора приближение методом РК четвертого порядка вызывает меньшие колебания скоростей подач (0.0276%) [36] и не требует расчёта второй производной функции С (и). Наиболее распространённый четырехстадийный метод РК будет рассчитываться на основании уравнения

и]+1 = и] +1 (кг + 2к2 + 2к3 + к4)к, (20)

где Л = Лз] - величина шага, к1, к2, к3, к4 - вспомогательные коэффициенты, рассчитываемые по формулам

ki = f(uj) = fi^

1 1 (21)

\\C'(uj)\\ ^x2+y2+ z2'

1 (22)

k2 = f(uj + 0.5k1Asj) = -T7-;

k3 = f(u,j + 0.5k2Asj) = t|—-

C'(uj + 0.5k1Asj)

1 (23)

\\C'(uj + 0.5 k2Asj)\\'

1 (24)

k4=f(Uj+k3ASj) = \Fiuj + k3Asj)\\.

Однако как показывают последние исследования [14,36], наибольшей точностью приближения параметра сплайна обладает метод РК второго порядка с компенсационной схемой приближения, где отклонение подачи составило всего 10-6 %. Связь последующего параметра сплайна uj+1 с текущим значением для данного метода записывается в виде уравнения

uj+1 = uj+1,temp + Auj+i, (25)

где uj+1temp - временный параметр, рассчитываемый по двухстадийному методу Рунге-Куты, который характеризуется уравнением

Asj (26)

uj+1,temp = uj + + k2).

Соответствующие коэффициенты из (26) определяются по формулам k = 1 k = 1 (27)

1 \\C'(uj)\y 2 \\C'(Uj + kiASj)W

При расчёте компенсационного параметра сплайна Auj+1 в идеальном случае должно выполняться равенство (17), т.е.

\\C(uj+i,temp + Auj+i) - C(uj)\\ = Asj. (28)

Используя разложение в ряд Тейлора первого порядка и свойство (17), запишем выражение для

расчёта uj+i, temp:

uj+1 = uj+1,temp + cTti Т (C(uj+1,temp + Auj+1) — C(.uj+1,temp}^>

С (uj+1, temp) (29)

С(и]+1£етр + Аи]+1) = С(и]+1£етр) + С (и]+1Хетр)Аи]+1

Подставляя (29) в (28) и возведя обе части в квадрат, получим следующее равенство:

||С(М]'+1, 1етр) + С'(и]+1,1етр)Аи]+1 - С(и])\\ = А$2. (30)

Сформируем квадратное уравнение вида

ААи2+1 + ВАи]+1 + 0 = 0, (31)

где коэффициенты А, В, О определяются по формулам А = \\С (и]+1{етр^ \\ ,

(32)

В = 2 (с'(и]+1,1:етр), {С(и]+1Хетр) — С(и]))),

О = ЦС(и]+1,сетр) - С(и]-)Ц -Аз].

Действительные корни уравнения (31) Аи]+11, Аи]+12 определяются по формуле

_-В± VВ2 - 4АО (33)

Аи]+112 = Та .

Учитывая, что метод обладает третьим порядком точности, можно принять значение О « 0. Подставляя его в (33), можно обнаружить, что

Аи]+12 « 0,а Аи]+12 « —В/А. Для обеспечения стабильности интерполируемых данных в (33) всегда используется относительно малое значение параметра

Аи]+1, т.е. первый корень. В случае, если корни мнимые, то значение компенсационного параметра принимается равной нулю. Таким образом, получим окончательное соотношение для Аи]+1 , которое затем подставляется в (25), в виде

-В + VВ2 - 4AD „ (34)

Аиц.Л =

4+1 = { 2A

0, В2 - 4AD < 0

, В2 - 4AD > 0

Можно заметить, что метод не требует определения второй производной, а первая производная рассчитывается три раза. Алгоритм демонстрирует в среднем трёхкратное снижение быстродействия по сравнению с методом разложения в ряды Тейлора второго порядка.

Хорошей точностью приближения также обладают методы поиска параметра сплайна с обратными связями, например, предсказывающе-исправляющий метод [31] (predictor-corrector), реализующий на первом этапе разностную схему аппроксимации производных, а на втором -схему коррекции параметра сплайна до того момента, когда ошибка не будет удовлетворять заданному допуску. Схема предиктор-корректора в некоторых случаях может выполняться быстрее, чем все описанные ранее в статье методы, однако время расчёта может сильно зависеть от подобранных параметров алгоритма и геометрии кривой, т.е. алгоритм не является численно стабильным. Данные недостатки были решены в алгоритме слежения за хордой (chord-tracking

algorithm) [30], в котором сначала применяется разложение в ряд Тейлора со схемой компенсации длины дуги, а затем выполняется уточнение скорости подачи при помощи итерационного алгоритма Ньютона-Рафсона. Метод слежения за хордой имеет квадратичную сходимость и уже после одной итерации максимальное отклонение скорости подачи имеет порядок 10-6 %, тогда как предиктор-корректор имеет значение около 0.1%.

3 Результаты вычислительных экспериментов

Для моделирования обработки использовался пространственный контур бабочки (рис. 4а, б), состоящий из 127 линейных сегментов, сформированных в CAM-системе 'TеММа-3D". В результате применения геометрического модуля сглаживания со вставкой кубических B-сплайнов исходный контур преобразовывался в линейно-сплайновую траекторию. Моделирование интерполятора осуществлялось в программном обеспечении, написанном на языке C++ в интегрированной среде разработки QtCreator. Численные эксперименты были посвящены исследованию линейно-параметрической интерполяции, где на сплайновых участках рассматривалось четыре наиболее распространённых метода прямого приближения параметра сплайна: аппроксимация рядами Тейлора второго порядка, аналогичного метода со схемой компенсации длины дуги кривой, метода РК четвертого порядка и двухстадийного метода РК с компенсационной схемой приближения. В экспериментах также учитывался алгоритм s-образного планирования подач с контролем рывка со следующими кинематическими ограничениями:

заданная скорость подачи = 166,667 мм/с (10 м/мин), центростремительное ускорение

2 2 «-» Аптах = 498 мм/с , тангенциальное ускорение АШах = 498 м/с , центростремительный рывок

3 «-» з о

1птах = 2000 мм/с и тангенциальный рывок ]£тах = 2000 мм/с . При построении B-сплайна

задавались следующие геометрические ограничения: максимально допустимая хордовая ошибка

8 = 0,005 мм, допуск ошибки аппроксимации сплайна £ = 0,1 мм, отношение длин сглаживания

с = 0,25.

По результатам экспериментов с учётом заданного такта управления Т3 = 400 мкс интерполируемый контур составлен из 56 864 точек, а продолжительность отработки траектории составляет «22с. (рис. 5а). Хордовая ошибка построения сплайна рассчитывалась по соотношению (3) и его максимальное значение составило 8тах ~ 6.1 • 10-6 мм (рис. 5б), а максимальная кривизна ктах « 206,1 мм-1 (рис. 5в). Как можно заметить, реальная ошибка 8тах оказалась намного меньше, чем 8, из-за преобладания других ограничений в скоростях подач, связанных с нормальными составляющими рывка и ускорения.

» -га о а »

Проекция координат X« ни

(б)

Рис. 4 Интерполируемая траектория тестового контура, проекция точек на плоскость ХУ (а),

проекция точек на плоскость Х2 (б).

(в)

Рис. 5 Результат применения алгоритма планирования подач для всей траектории (а), зависимость кривизны траектории к (б), хордовая ошибка построения сплайна 8 (в).

Первый проход с аппроксимацией параметра сплайна рядами Тейлора втрого порядка демонстрирует максимальные колебания скоростей подач в пределах «0.32% (рис. 6а). Практически такие же значения показывает алгоритм с компенсацией длины дуги, где разница между графиками видна только при увеличении масштаба (рис. 6б). Наибольшая амплитуда отклонений наблюдается в местах перехода сплайнового и линейного участков, а метод с компенсацией длины дуги имеет лучшую точность при большей кривизне сплайна.

Рис. 6 Сравнение колебаний скоростей подач при использовании рядов Тейлора второго порядка (зеленый) с аналогичным методом при компенсации длины дуги (синий) (а), при

увеличении отдельной области (б).

Использование метода РК четвертого порядка показывает отклонения в скоростях подач до 0.036% (рис. 7а). При этом значительная ошибка наблюдается в зоне наибольшей кривизны. Резкие колебания сохраняются и при сопряжении линейного и сплайнового участков (рис. 7б), хотя ошибка не превышает « 6,5 • 10"5 %. Последний четвертый проход контура, в котором используется метод РК с компенсационной схемой приближения, демонстрирует наилучшие

у

результаты, где флуктуации равны « 1,681 • 10" % (рис. 7в). Как и в случае с разложением в ряд Тейлора, максимальные колебания встречаются в линейно-сплайнового переходах, а к середине кривой эти колебания уменьшаются (рис. 7г).

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

(г)

Рис. 7 Колебания скоростей подач при использовании алгоритма РК четвертого порядка (а) и при увеличении отдельной области (б), метод РК с компенсационной схемой приближения

(в), при увеличении отдельной области (г).

4 Заключение

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

Анализ наиболее распространённых методов параметрической интерполяции, используемых в системах управления станков и промышленных роботов, а также сравнение результатов экспериментов показывают, что качество отработки контура для одной и той же кривой зависит от применяемого алгоритма параметрической интерполяции. Наименьшие колебания скорости подачи наблюдаются при использовании алгоритма РК второго порядка с компенсационной схемой приближения, которая полностью удовлетворяет требованиям к высокоскоростной обработке и рекомендуется для использования в системах управления станков и роботов.

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

Благодарности

Исследование выполнено за счёт гранта Российского научного фонда N0 21-79-10392, https://rscf.ru/project/21 -79-10392/.

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

[1] Brecher С., Lange S., Merz M., Niehaus F., Wenzel C., Winterschladen M. NURBS based ultra-precision free-form machining. CIRP Annals, vol. 55 (1), pp. 547-550, 2006.

[2] U.S. Product Data Association. Initial Graphics Exchange Specification, IGES 5.3, 1996.

[3] Pratt. M. J. Introduction to ISO 10303-the STEP standard for product data exchange. Journal of Computing and Information Science in Engineering, vol. 1, no. 1, pp. 102-103, 2001.

[4] Zhong W.B., Luo X.C., Chang W.L., Cai Y.K., Ding F., Liu H.T., Sun Y.Z. Toolpath Interpolation and Smoothing for Computer Numerical Control Machining of Freeform Surfaces: A Review. International Journal of Automation and Computing, vol. 17(1), pp.1-16, February 2020.

[5] Абдуллин Т.Х., Харьков М.А. Алгоритм опережающего просмотра для системы ЧПУ с применением трапецеидальных законов разгона/торможения. "Машиноведение и инновации. Конференция молодых учёных и студентов", 2017 6-8 декабря.

[6] Zhong W.B., Luo X., Chang W., Ding F., Cai Y. A real-time interpolator for parametric curves. International Journal of Machine Tools and Manufacture, vol. 125, pp.133-145, 2018.

[7] Ernesto C., Farouki R. T. A. High-speed cornering by CNC machines under prescribed bounds on axis accelerations and toolpath contour error The International Journal of Advanced Manufacturing Technology, vol. 58, pp. 327-338, 2012.

[8] Yeh S-S., Hsu P-L. Adaptive-feedrate interpolation for parametric curves with a confined chord error. Computer-Aided Design, vol. 34(3), pp. 229-237, 2002.

[9] Nam S-H., Yang M-Y. A study on a generalized parametric interpolator with real-time jerk-limited acceleration. Computer-Aided Design, vol. 36, pp. 27-36, 2004.

[10] Zhao H., Zhu LM., Ding H. A real-time look-ahead interpolation methodology with curvature-continuous B-spline transition scheme for CNC machining of short line segments. International Journal of Machine Tools & Manufacture, vol. 65, pp. 88-98,2013.

[11] Beudaert X., Lavernhe S., Tournier C. 5-axis local corner rounding of linear tool path discontinuities, International Journal of Machine Tools & Manufacture, vol. 73, pp. 9-16, 2013.

[12] Заруднев А.С., Илюхин Ю.В. Повышение производительности мехатронных комплексов лазерной обработки на основе зависимости контурной погрешности от параметров движения и исполнительных систем. Известия Самарского научного центра Российской академии наук, Том 9, №3 (21), с. 758 - 764, 2007.

[13] Заруднев А.С., Илюхин Ю.В. Повышение производительности лазерных комплексов на основе прогноза контурной ошибки. Изд-во «Новые технологии», М.: Мехатроника, автоматизация, управление, № 9, с.52-56, 2010.

[14] Huang J., Du X., Zhu L-M. Real-time local smoothing for five-axis linear toolpath considering smoothing error constraints. International Journal of Machine Tools and Manufacture, vol. 124, 67-79, 2018.

[15] Bi Q-Z., Wang Y-H., Zhu L-M., Ding H. A Practical Continuous-Curvature Bezier Transition Algorithm for High-Speed Machining of Linear Tool Path. Intelligent Robotics and Applications: 4th International Conference, ICIRA, Aachen, Germany, December 6-8, Proceedings, Part II, pp. 465-476, 2011.

[16] Fan W., Lee C-H., Chen J-H. A real-time curvature-smooth interpolation scheme and motion planning for CNC machining of short line segments. International Journal of Machine Tools and Manufacture, vol. 96, pp. 27-46, 2015.

[17] Farouki R.T. Construction of G2 rounded corners with Pythagorean-hodograph curves. Computer Aided Geometric Design, vol. 31, pp. 127-139, 2014.

[18] Wang J.B., Yau H.T. Real-time NURBS interpolator: application to short linear segments. The International Journal of Advanced Manufacturing Technology, vol. 41, pp. 1169-1185, 2009.

[19] SIEMENS. SINUMERIK 840D sl/828D Basic Functions Function Manual. [Электронный ресурс]. - Режим доступа: https://cache.industry.siemens.com/dl/files/431/109476431/att 844512/v1/FB1sl 0115 en , свободный. - (дата обращения: 30.06.2020).

[20] FANUC Series 30i/31i/32i/35i-MODEL B [Электронный ресурс]. - Режим доступа: https://www.fanuc.co.jp/en/product/catalog/pdf/cnc/FS30i-B(E)-09a.pdf, свободный. - (дата обращения: 30.06.2020).

[21] Chen Z.C., Khan M.A. Piecewise B-Spline Tool Paths With the Arc-Length Parameter and Their Application on High Feed, Accurate CNC Milling of Free-Form Profiies. Journal of Manufacturing Science and Engineering, vol. 134 (3), 2012.

[22] Основы кривых Безье [Электронный ресурс]. - Режим доступа: https://pomax.github.io/bezierinfo/ru-RU/index.html , свободный. - (дата обращения: 15.06.2020).

[23] Зеленский, А.А., Абдуллин Т.Х., Илюхин Ю.В., Харьков М.А. Высокопроизводительная цифровая система на основе ПЛИС для управления движением многокоординатных станков и промышленных роботов. "СТИН" (СТанки Инструмент), №8, c. 5-8, 2019.

[24] Зеленский, А.А. Харьков М.А., Ивановский С.П., Абдуллин Т.Х. Высокопроизводительная система числового программного управления на базе программируемых логических интегральных схем. Вестник Воронежского государственного технического университета, Том 14, №5, c. 8-12, 2018.

[25] Зеленский, А.А., Шадрин Н.Г., Абдуллин Т.Х., Харьков М.А. Высокоскоростная промышленная сеть реального времени киберфизических систем. Вестник компьютерных и информационных технологий, №11, c. 46-52. 2019.

[26] Yuen A. Spline interpolation and contour error pre-compensation for 5-axis machining Dissertation. [Электронный ресурс]. - Режим доступа: https://open.library.ubc.ca/soa/cIRcle/collections/ubctheses/24/items/1.0073966 , свободный. - (дата обращения: 20.08.2020).

[27] Koren Y., Lo C. C., Shpitalni M. CNC interpolators: Algorithms and analysis. Manufacturing Science and Engineering. Manufacturing Science and Engineering, vol. 64, pp. 83-92, 1993.

[28] Bhattacharjee. B., Azeem A., Ali S., Paul S. Development of a CNC interpolation scheme for CNC controller based on Runge-Kutta method. International Journal Computer Aided Engineering and Technology, vol. 4, № 5, pp. 445-464, 2012.

[29] Cheng M-Y., Cheng M-Y., Tsai M-C., Kuob J-C. Real-time NURBS command generators for CNC servo controllers. International Journal of Machine Tools and Manufacture, vol. 42, pp. 801-813, 2002.

[30] Zhao H., Zhu L., Ding H. A parametric interpolator with minimal feed fluctuation for CNC machine tools using arc-length compensation and feedback correction. International Journal of Machine Tools & Manufacture, vol. 75, pp. 1-8, 2013.

[31] Tsai M.C., Cheng C.W. A real-time predictor-corrector interpolator for CNC machining. Transactions of ASME Journal of Manufacturing Science and Engineering, vol. 125, pp. 449460, 2003.

[32] Heng M., Erkorkmaz K. Design of a NURBS interpolator with minimal feed fluctuation and continuous feed modulation capability. International Journal of Machine Tools and Manufacture, vol. 50, pp. 281-293, 2011.

[33] Zhang L.B., You Y.P., He J., Yang X. F. The transition algorithm based on parametric spline curve for high-speed machining of continuous short line segments. The International Journal of Advanced Manufacturing Technology, vol. 52, pp. 245-254, 2011.

[34] Зеленский, А.А., Абдуллин Т.Х., Алепко А.В. Особенности построения в реальном времени s-образной кривой разгона/торможения при кусочно-линейной интерполяции поверхностей сложной формы. Робототехника и техническая кибернетика, том 9, № 3, с. 186-195, 2021.

[35] Farouki R.T., Sakkalis T. Real rational curves are not unit speed. Computer Aided Geometric Design, vol. 8, pp. 151-157, 1991.

[36] Jia Z-Y., Song D-N., Ma J-W., Hu G-Q., Su W-W. A NURBS interpolator with constant speed at feedrate-sensitive regions under drive and contour-error constraints. International Journal of Machine Tools and Manufacture, vol. 116, pp. 1-17, 2017.

Research of parametric interpolation approaches for smoothing a

piecewise linear toolpath

Zelensky A.A.1' , Abdullin T.H.1, , Dubovskov V.V.1, , Ilyukhin Yu.V.1, v

1mMSTU "STANKIN"

*zelenskyaa@gmail.com **everestultimate@yandex.ru ***dtntll2l@gmail.com ****ilyv-2@mail.ru

Abstract. The paper considers the basic methods of parametric interpolation of tool trajectory during contour machining of products having complex geometric surface shape. The paper applies an approach based on inserting cubic B-splines into the initial piecewise linear trajectory and implementing linear and parametric interpolation algorithms which allow increasing machining productivity. When there is no explicit analytical dependence spline parameters on its length of the corresponding trajectory segment parametric interpolation is performed. Therefore, numerical methods are used for formation of interpolated points along the curve, which have different accuracy of approximation of the spline parameter and have a significant impact on the quality of machining of the product due to the unstable feed rate. A number of numerical experiments have been conducted to select the most effective parametric interpolation algorithm for a trajectory that has a geometric continuity of G2. The best results are shown by the second-order Runge-Kutta algorithm with a compensating approximation scheme, which is the best choice for control systems of machine tools and industrial robots for high-precision and high-speed machining of parts with complex geometries.

Keywords: parametric interpolation, Runge-Kutta algorithm, Taylor series expansion, trajectory smoothing, motion controller, high-speed processing.

Acknowledgements This work is supported by the Russian Science Foundation under grant No 2l-79-10392, https://rscf.ru/project/2l-79-l0392/