CC BY
90
Н.К.Кузнецов, А.Ю.Перелыгина
Разработка программного обеспечения задач динамики и управления движением манипуляционных роботов
Проблема динамического синтеза управляемых машин с упругими звеньями связана с необходимостью составления уравнений движения и решения прямых и обратных задач кинематики и динамики. Решение этой проблемы в манипуляционных роботах осложняется большим количеством степеней подвижности, переменностью параметров и кинематической избыточностью исполнительных механизмов. Для ее решения были разработаны автоматизированные процедуры и создано большое число программных комплексов, позволяющих составлять дифференциальные уравнения исполнительных механизмов роботов и решать прямые и обратные задачи кинематики и динамики [1,2],
Однако известные комплексы были реализованы на устаревшем программном обеспечении типа MS-DOS (версии 3-5), MacOS и таких файловых менеджеров, как Norton Commander, Volkov Commander, Far Manager, имеющих ограниченный графический интерфейс, малую гибкость и адаптивность к требованиям пользователей, и вследствие, этого плохую наглядность и сложность в выводе и обработке полученной информации. Указанные недостатки ограничивали возможности их использования для решения задач управления движением манипуляционных роботов с учетом упругих свойств конструкции в режиме реального времени. В настоящей статье описывается программный комплекс для решения прямых и обратных задач кинематики и динамики и управления движением манипуляционных роботов, позволяющий учитывать упругость звеньев и определять потребные управляющие силы и моменты [3].
Для решения прямой задачи кинематики декартовые координаты, линейные скорости и ускорения точек цепи в разложении на заданные оси в разработанном комплексе программ определяются по следующим зависимостям:
xj=xj(xnp,q); rj=rj(xnp,q,q); (D
где Xj - вектор декартовых координат точки х. в системе /; xi - заданный вектор декартовых координат точки цепи
в системе /'; V] - вектор линейной скорости точки X, в системе /; p,q,q,q - векторы параметров, обобщенных
координат, скоростей, ускорений.
При решении прямой задачи используется последовательное вычисление координат точки xt во всех системах с номерами г < j, пока текущий индекс г Ф j. Рекуррентные формулы при этом имеют вид
(2) (3)
X;
(г-1)
мх,
для случая поворота систем /', /-1,
X:
(г-1)
= X; + L
г-1 ?
для случая переноса.
В (3) - вектор, содержащий одну ненулевую компоненту; перенос х: Ьг_{ ,0,0); перенос у:
Ьг1 - (0,/гЧ,0) ; перенос т. Ьг_л — (0,0,; = цг.А - если перенос переменный; = /г_1 - если постоянный.
Матрицы Мг в (2) имеют вид
( 1 0 0 cos /,. 0 sin N ( 0 0 °1
0 cos sin 7 5 0 1 0 5 cos f, - sin / 0
sin /, cos J ч - sin / ; о cos /,у f, cos fi h
соответственно осям X, у, г относительно поворота систем - /-1; fi = Я, - если поворот переменный; f:
если постоянный. Выбор (2), (3) и величины f 1 выполняется в зависимости от вида перехода между системами координат. Выражения для скоростей и ускорений можно получить дифференцированием (2), (3) по времени
У1Г'Х) =Цх1г) +Мгх(р;
-Мгх{р +2Му\г)+МгУ\г\
для случая поворота,
у С-1) -
г
К(г) + К;
-- + IV
(4)
(5)
г '
для случая переноса.
Алгоритмы решения обратных задач кинематики о положении строятся на основе минимизационных методов. В отличие от рассмотренных выше прямых задач кинематики решение обратных задач может быть неоднозначным. Так, если число условий т, наложенных на положение системы, связанной со схватом, больше числа обобщенных координат п, решений может не быть совсем; если п - т, то число решений, если они есть, будет ограничено; для случая
п > т имеем множество решений, принадлежащих некоторой области О. Следовательно, при решении обратной задачи кинематики существуют следующие ограничения:
1) на обобщенные координаты механизма не наложено никаких ограничений;
2) обобщенные координаты должны удовлетворять ограничениям вида цн <д<дн\
3) в инерциальной системе координат задана только точка хс[а, Ь, с), с которой требуется совместить начало системы, связанной со схватом;
4) дополнительно к координатам точки заданы ещё и углы ориентации, связанной со схватом системы осей
(гЛг)-
При составлении дифференциальных уравнений движения исполнительные механизмы представляются в виде звеньев, соединенных между собой кинематическими парами или сосредоточенными упругостями, причем сами звенья могут состоять из одного и более тел, связанных упругими связями, а последние интерпретируются как кинематические пары, в которых действуют упругие и демпфирующие силы. Такой прием позволяет с достаточной точностью представлять упругие звенья исполнительных механизмов, для чего они разбиваются на требуемое число твердых тел, соединенных упругими элементами, движения в которых есть деформации изгиба, кручения и растяжения-сжатия. Динамика манипулятора описывается системой обыкновенных дифференциальных уравнений
Ои+йг+ОсХя)
где А(д) - симметричная, положительно определенная матрица коэффициентов инерции (п х т); В(д)д - вектор скоростных сил; (¿и ,0.? ->(1(;{д) - векторы сил управления, внешних сил и сил тяжести, приведенные к обобщенным координатам; — с Ад - упругие силы, приведенные к осям упругих деформаций; - к Ад - силы демпфирования, приведенные к осям упругих деформаций; - активные силы, приведенные к осям упругих деформаций; Сг, (д) - силы тяжести, приведенные к осям упругих деформаций.
Для вычисления элементов матрицы инерции и вектора демпфирующих сил используется алгоритм, основанный на принципе Гаусса и уравнениях Аппеля, согласно которому левую часть системы (6) можно получить, определяя частные производные дС!дд, где в - энергия ускорения цепи [1]. Так как главные центральные оси инерции звеньев параллельны осям связанных систем координат, то энергия ускорения записывается следующим образом
0=X +X X + ХМ(7)
(6)
Тх 2
М №
ы
При этом предусмотрены два режима получения дифференциальных уравнений движения - с учетом и без учета упругих деформаций, их скоростей и ускорений. Интегрирование уравнений движения производится методом Рунге-Кутта четвертого порядка с автоматическим выбором шага интегрирования по заданной погрешности, Для этого система уравнений (6) приводится к нормальной форме Коши
I Р = А-] (Я)[йи +(2Р+ ес (д) - В(д, р)1,
где p = q - вектор обобщенных скоростей; A l(q)- матрица, обратная к А, в уравнении (6), для чего можно использовать более простую и быструю, чем обращение матрицы, процедуру решения (6) как системы линейных алгебраических уравнений относительно неизвестных ускорений q ,
Для решения обратных задач динамики, заключающихся в определении необходимых управляющих сил по заданным законам изменения координат, движение исполнительного механизма робота разделяется на программное и колебательное. При этом управляющие силы, реализующие программное движение, находятся с помощью уравнений (6) по заданному уравнению траектории, без учета влияния колебательного движения
Управляющие силы, реализующие программное движение, находятся с помощью уравнения (6) по заданному уравнению траектории, без учета влияния колебательного движения
Ql = )q* + B(q* )q* - gF - QG (q* ), (9)
где - вектор координат программного движения; Qu - вектор программных управляющих сил.
Дополнительные силы, осуществляющие компенсацию упругих колебаний, определяются из уравнения колебательного движения
A(q* )A'q + В Ас/ + cAq = - A{q* )'q*. (10)
Здесь A q - вектор упругих деформаций; В - диагональная матрица коэффициентов демпфирования; С - диагональная матрица коэффициентов жесткости [4]. Вектор ускорений программного движения <jf* в уравнении (10) находится на основе уравнений (6) с учетом взаимовлияния программного и колебательного движений
Ч* = A'1 (q*)[Q„ + AQU +Qf+ Qg(q, ) - B(q, )]qt - BAq - cAq, (11)
где AQu - дополнительные силы, компенсирующие упругие колебания.
Программный комплекс организован на алгоритмическом языке Фортран с компилятором Microsoft Fortran Power-Station 4,0 и имеет модульную структуру [5]. Он основан на использовании как стандартных программ, входящих в библиотеки математического обеспечения ЭВМ, так и специальных программ, учитывающие специфические особенности динамики исполнительных механизмов роботов. Все геометрические и массоинерционные характеристики, а также вид расчетной схемы задаются пользователем в виде диалога. С целью эффективной организации счета и удобства пользования комплексом программ данные, которые изменяются в процессе счета или задаются пользователем, отделены от данных, которые остаются неизменными для моделируемой кинематической цепи, Первую группу данных названы данными задач, вторую - данными моделей.
Кинематическая цепь исполнительного механизма описывается следующим массивом параметров
М\ = (Ск>Сп>С*>Ст>Ср>К>к*>Ь><*>(1н>Яв)' где ск - число связанных систем координат (цифр в коде цепи); сп - число обобщенных координат основного движения (цифр кода 1 <кц <6); с5 - число обобщенных координат упругих деформаций (цифр кода 11 ¿кц £16); ст - число масс (цифр кода <-10); с - число параметров (цифр кода <0); кц - код кинематической цепи; ks - код зависимости обобщенных координат (с число элементов сп ); b,d - коэффициенты связи; qH,qB - верхние и нижние пределы изменения обобщенных координат.
При решении прямой задачи кинематики пользователем задаются характеристики относительного движения такие, как обобщенные координаты, обобщенные скорости и ускорения. Результатом работы комплекса является определение положения, ориентации, линейной и угловой скоростей и ускорений звеньев исполнительного механизма робота в инерционной системе координат. И, наоборот, при решении обратной задачи кинематики исходными данными являются абсолютные координаты положения схвата, линейные скорости и ускорения звеньев или отдельного звена механизма. В результате определяются характеристики относительного движения, соответствующие заданному положению и ориентации схвата.
Для решения прямых задач кинематики предназначены подпрограммы KINX, KINV, KINW. Подпрограмма KINX служит для определения декартовых координат точки в системе с номером N, по её координатам в системе с номером М; её вызов имеет вид CALL K!NX(Q, SK, DX, M, N, IPZ). Здесь используются следующие входные данные: M - номер системы, в которой заданы координаты точки; N - номер системы, в которой требуется вычислить её координаты; SK(3) - массив координат (X,y,z) точки в системе M; IPZ - ключ, если IPZ=0, вычисления выполняются без учета упругих деформаций, иначе - с учетом; Q(NSUM) - массив обобщенных координат; Требуемые общие области модели состоят из /KODC/, /VIKP/, /PARIM/. Выходными параметрами являются DX(3) - массив координат точки в системе N. Значения M и N могут быть любыми 8 пределах от 0 до KS. Для формирования вектора индексов по коду цепи предназначена под-
программа KOD, без параметров, Её входные данные содержатся в общих областях /CONSM/ и /KODC/; результат помещается в массив IKP блока /VIKP/'.
Подпрограмма KINV выполняет вычисления декартовой и линейной скорости заданной точки цепи. Её вызов имеет вид CALL KINV(Q, V, SK, DX, DV, M, N, IPZ). Кроме данных, уже описанных в вызове KINX, требуется задать V(NSUM) -вектор обобщенных скоростей цепи. Выходными параметрами являются: DX(3) - массив декартовых координат точки в системе N; DV(3) - массив линейной скорости точки относительно системы N в разложении на её оси, Предполагается, что заданная точка неподвижна относительно системы М, кроме того, в этом случае требуется М > N, в противном случае вычисления не выполняются,
Для вычисления линейного ускорения точки предназначена подпрограмма KINW, имеющая вызов CALL KINW(Q, V, W, SK, DX, DV, DM, N, M, IPZ). Здесь в дополнение к данным KINV требуется задать W(NSUM) - массив обобщенных ускорений. Результат вычисления линейного ускорения заданной точки в её движении относительно системы N в разложении на оси х, у, z этой системы помещается в массив DW(3); при этом, как и в KINV, требуется М > N. Подпрограммы KINX, KINV, KINW осуществляют последовательные вычисления координат, скоростей и ускорений в системах М - 1, М - 2 и т.д., пока текущий индекс системы не будет #!.
Для удобства в разработанном программном комплексе программы KINX, KINV, KINW объединены в три подпрограммы для поворотов относительно х, у, г соответственно в KINX, KINY, KINZ. Вызов последних имеет вид
CALL KINX(Q0, V0, АО, Ql.Vl, AI, t);
CALL KINY(Q0, VO, AO, Ql.Vl, AI, t);
CALL KINZ(Q0, VO, AO, Ql.Vl, AI, t), где QO, Ql— массив начальных и полученных значений координат в разложении на оси соответственно; VO, VI - массив начальных и полученных значений скоростей в разложении на оси координат соответственно; АО, AI - массив начальных и полученных ускорений в разложении на оси координат соответственно; t - текущий момент времени.
Комплекс программ решения обратных задач кинематики о положении включает в себя подпрограммы: BACK - головная подпрограмма; KINOX, KINOY, KINOZ - те же подпрограммы, что и KINX, KINY, KINZ , только для определения координат методом подбора обобщенных координат. Вызов подпрограмм имеет вид
CALL BACK(KP,FP,M, IER),
Для решения прямой задачи динамики в качестве исходных данных используются характеристики привода, обобщенные координаты, скорости, и ускорения, условия интегрирования (начальный момент временного интервала, конечный момент интервала интегрирования, шаг интегрирования, допустимая погрешность на шаге интегрирования). В случае решения обратной задачи динамики исходными данными служат конечные или дифференциальные уравнения изменения упругих координат.
Подпрограмма APPEL вычисления матрицы А и вектора В уравнений динамики (6) для механизмов с произвольной кинематикой основана на следующем алгоритме.
■ Производится вычисление энергии ускорений G0 звеньев кинематической цепи при нулевом векторе обобщенных ускорений q = 0. Одновременно находятся значения абсолютных координат и скоростей центров масс звеньев, которые записываются в рабочий массив WR, имеющий размер
1 = 4^(5, + /,) + «,
где г - число масс; п - число обобщенных координат; •Sj, Ц - число переменных и постоянных поворотов в коде цепи,
считая от стойки до системы, в которой расположена /-я масса.
■ Вычисление диагональных элементов матрицы инерции и вектора демпфирующих сил сводится к вычислению энергии ускорений.
■ Для определения недиагональных коэффициентов а у вычисление G выполняется для вектора qx, = (0,...1,...0Д,...0) , где единицы стоят на /-м и /-м местах при нулевых значениях обобщенных скоростей.
Подпрограмма APPEL может работать в двух режимах: в первом вычисляется матрица инерции и вектор демпфирующих сил программного движения, при этом обобщенные координаты и скорости упругих деформаций игнорируются; во втором режиме - полная матрица инерции и вектор демпфирующих сил с учетом всех обобщенных координат цепи.
Для вычисления вектора обобщенных сил предназначен модуль FORCE, вызов которого имеет вид
CALL FORCE (A,B,Qc,t).
Здесь входными данными являются матрица инерционных характеристик А и вектор сил В , а выходом - вектор обобщенных сил Qc.
Процедура интегрирования реализуется подпрограммой RUNGE, которая вызывает подпрограммы ACCEL приведения (5) к нормальной форме (8) и подпрограмму PRINM вывода результатов интегрирования на принтер. Подпрограмма ACCEL вызывает подпрограммы MODELM и MODELQ моделирования исполнительного механизма и силовой нагрузки. Подпрограммы PRTNM, MODELQ, MODELM составляются пользователем. Для удобства пользователя вызов подпрограммы RUNGE имеет следующий вид
RUNGE (fq.vq.t.tn.tk.a.b).
Здесь fq - массив обобщенных координат; vq - массив обобщенных скоростей; wq - массив обобщенных ускорений; I, in, tk - текущий, начальный и конечный моменты времени соответственно; a, b - коэффициенты динамической модели. Затем следует:
■ вычисление обобщенных сил управления (в функции времени Г обобщенных координат q, скоростей V);
■ вычисление сил тяжести и приведение их к обобщенным координатам;
■ вычисление обобщенных сил от активного внешнего нагружения; суммирование результата и запись его в массив QS,
Для решения обратной задачи динамики, заключающейся в определении необходимых управляющих сил по заданным законам изменения упругих координат, используется следующий алгоритм:
■ по заданному уравнению траектории вычисляются значения обобщенных координат q*(t), скоростей </* (0 и ускорений <7* (0 программных движений;
■ одновременно на основе заданных зависимостей осуществляется вычисление упругих координат Дq(f) и их производных;
■ по известным q *, q * находятся значения элементов матрицы инерции A(q*), вектора демпфирующих сил B(q* )qt и сил тяжести в (9);
■ определяются программные управляющие силыОц, для чего найденные значения q(t\q(t\A(c[*) и
B{q*) подставляются в (9);
■ значения программных управляющих сил, упругих координат и их производных вводятся в (11);
■ с помощью уравнений (11) и (10) определяются силыЛ»2ц1 осуществляющие компенсацию упругих колебаний;
« осуществляется переход к следующему моменту времени t¡+i=t¡+At, где At - шаг интегрирования;
■ вычисления продолжаются до истечения заданного интервала времени [0,Т],
Продемонстрируем возможности использования предложенного комплекса программ для анализа свободных колебаний исполнительного механизма робота с пятью степенями подвижности. Расчетная схема этого манипулятора приведена на рис.1.
Упругие кинематические пары на этом рис. помечены знаком £ Упругость звена ^ приведена в точку 1 и представлена в виде упругого сферического шарнира, звена 12 - в точку 2, в которой выступает в виде упругой сферической пары с пальцем. Сами звенья ^ и 12 при этом считаем абсолютно жесткими. Этот механизм имеет 10 степеней свободы; его конфигурация определяется обобщенными координатами с^ + с/5; относительно положения равновесия происходят малые колебания (деформации) в упругих шарнирах 1 и 2, описываемые обобщенными координатами * Ац5.
Вывод результатов работы программного комплекса производится в виде отдельного окна или графиков, На рис. 2 представлено окно решения задачи, в котором вычисляются матрица инерции программного движения, полная матрица инерции; матрица инерции колебательного движения и векторы скоростных сил; спектр частот собственных колебаний; матрица коэффициентов собственных форм; относительные амплитуды и начальные фазы,
ш modd Ш: 11ДШИ0Ш
\яёшияюШ 1 ш ш |И-И|вЖ1
МАТРИЦА ИНЕРЦИИ ОСНОВНОГО ДВИЖЕНИЯ
4.239653 0 . 0. 0.194709 463119 0.0 3.896942 2.408742 0. 0 -0.264591 0.000001 2.408742 1. 583666 0.0 -0.187416 0.194709 0.0 0.0 0.5 0.0 0.463119-0.264591 -0.187416 0.0 0.13«
ПОЛНАЯ МАТРИЦА ИНЕРЦИИ 4.235798 -0.011690 -0.007226 0.195089 0.463452 -0.011690 3 896940 2.438731 0.0 -0.264591 -0.007226 2.408731 1.583665 0.0 -0.187416 0.196089 0.0 0.0 0.5 0.0 0.463462 -0 264591 -0.187416 0.0 0.134800 -1.360496 0 022614 0.01 4001 0.459941 -0.222220 4.235783 -0.011650-0.0072260. 196089 0.463462 -0.001075 4.093945 2.571306 0.000587-0 278409 -0.645843 0.000002 0.0 0.459944 0.149746 -0.011690 3.896933 2.408736 0.0 -0.264590
МАТРИЦА ИНЕРЦИИ КОЛЕБАТЕЛЬНОГО ДВИЖЕНИЯ 2.236172-1.360496 0.020423 1.858138 0.022514 -1.360496 4.235798 -0.001075-0.645843 -0.011650 0.033423X1.001075 4.471931-0 001938 4 093945 1.858138-0.645845-0.001938 1.661143 0.С00002 0.022614-0.011650 4.093946 0.000002 3.89694]
ВЕКТОР СКОРОСТНЫХ СИЛ
0. 153436-0 271275 0.050416 0.0988090.019335
МАТРИЦА КОЭФФИЦИЕНТОВ СОБСТВЕННЫХ ФОРМ 1.0 10 1.0 1.0 1.0 0.159420 2.236942 0.777431 -1. 566849 -0.891978 -0.0С3727 0.003368 -227.5925 -0.006871 149.9349 -1 067851 2.407253 -1.123735 1. 488757 0.342313 -0.001414 XI. 003930 241.7648 -0.009190 403.2692 }ЕЯ=0 0 3 0 3
АМПЛИТУДЫ- 0.000651 0.003450 0. ООООСМ
0. 011438 0.000046 ФАЗЫ= -1. 234578 0.478752 2.1333244 -0.7752880.115420
Рис. 2. Результаты работы программного комплекса
Для графического представления результатов работы (рис. 3) использована стандартная программа EXCEL пакета MICROSOFT OFFICE, для чего в процессе вычисления результаты записываются в файл с расширением '.ехГ,
Рис. 3. Графическое представление результатов работы программного комплекса
Так как современный фортран обеспечивает совместимость с другими языками, возможен графический вывод результатов и с помощью других языков и программ, например, Visual С++, Visual Basic, Ассемблер, MathCAD и т.д.
Предложенный программный комплекс позволяет решать задачи динамики и управления движением манипуляцион-ных роботов с учетом упругих свойств звеньев в режиме реального времени. Этот комплекс может быть использован как в системах их автоматизированного расчета и проектирования, так и в системах программного управления. Его применение позволит повысить быстродействие и точность работы манипуляционных роботов различного назначения.
Библиографический список
1. Попов, Е. П. Манипуляционные роботы. Динамика и алгоритмы. / Е. П. Попов, А, Ф. Верещагин, С. Л, Зенкевич. - М. : Наука, 1978. - 400 с.
2. Бутырин, С. А. Пакет прикладных программ по моделированию и исследованию кинематики и динамики манипуляционных роботов / С, А, Бутырин [и др.] II Алгоритмы и программы : инф. бюлл, Госуд. фонда алгор, и прогр. СССР, - 1986, - № 3 (72). - С. 43.
3. Свидетельство об отраслевой регистрации разработки № 5641 код по ЕСПД ,03524577, 01299-01. Программный комплекс д\я исследования динамики и управления движением исполнительных механизмов роботов/ Н. К, Кузнецов, А. Ю, Перелыгина ; Иркут. техн. ун-т. - ГР № 50200600134 ; зарег, 16.02,06, - М. : ОФАП, 2005, - 5с.
4. Елисеев, С, В, Управление колебаниями роботов / С, В, Елисеев, Н. К. Кузнецов, А. В. Лукьянов. - Новосибирск : Наука. Сиб. отд-ние, 1990, - 320 с.
5. Бартеньев, О, В, Современный Фортран / О, В. Бартеньев. - М.: ДИАЛОГ- МИФИ, 2000, - 449 с.
Д.А.Журавлев, С.И.Ключников, Н.А.Карелина
Влияние интенсивности фрезерования на изменение пространственной формы маложестких авиационных панелей
Проблема обеспечения точности маложестких деталей летательных аппаратов при проявлении остаточных деформаций после механической обработки, относится к задачам от решения которых во многом зависят эксплуатационные показатели летательных аппаратов. Это связано с формированием в процессе механической обработки внутренних напряжений, влияющих на такие эксплуатационные показатели детали как: долговечность, надежность, статическая и усталостная прочность, релаксационная стабильность,
Одним из путей улучшения качественных характеристик деталей авиационного назначения является математическое прогнозирование и последующее технологическое управление пространственной формы маложестких деталей в процессе механической обработки.
Наиболее перспективным представляется создание пространственной модели деформирования тонкостенной детали на основе изменения технологических параметров фрезерования. При этом пространственное изменение конфигурации детали обусловлено воздействием различных сочетаний факторов технологического управления, влияющих на формирование различных по направлению пространственных тензоров деформаций.
Под технологическим управлением, в данном случае, понимаем регулирование комплексного параметра интенсивности фрезерования, включающего, в том числе, очередность фрезерования участков, на которые условно поделена деталь при ее дискретном рассмотрении и разгрузки заготовки от термических напряжений высверливанием несквозных отверстий.
Для оптимизации фрезерования маложестких деталей по критерию минимальных пространственных деформаций предложен следующий алгоритм:
1. Дискретное, поэлементное представление маложесткой детали.
2. Создание математического аппарата пространственных поэлементных деформаций как функции интенсивности фрезерования и разгрузки высверливанием несквозных отверстий,
3. Объединение отдельных элементов изделия в единую деталь с учетом функции взаимовлияния элементов и определение результирующей деформации.
Кратко рассмотрим каждый из перечисленных этапов.
Использование традиционных математических методов для расчета остаточных деформаций деталей летательных аппаратов не всегда позволяет получить удовлетворительные результаты, поэтому для определения поводок маложестких деталей предлагается использовать дискретный подход [1,2], который заключается в разбивке детали, на рад участков, при этом в дальнейшем, деформации каждого из участков (элементов) рассчитываются индивидуально.
