Текст работы

Математическое и имитационное моделирование, исследование устойчивости и стабилизация движения транспортных средств. Двухколесный велосипед

Совершенствование работающих и проектирование новых транспортных средств, сельскохозяйственной техники и других аналогичных механизмов является одной из основных тенденций развития современного машиностроения. Важнейшей задачей в этой области является построение математических моделей с целью исследования устойчивости и стабилизации движения этой техники.
В данной работе в качестве объекта исследования был выбран двухколесный велосипед. Из-за своей относительной простоты двухколесные транспортные средства являются более удобными для разработки математических моделей, исследования устойчивости и стабилизации движения. Обладая большей маневренностью по сравнению с трех и четырехколесными транспортными средствами, двухколесные, в том числе и роботы, представляют определенный интерес и с практической точки зрения.
Целью данного исследования является: 1) построение математической модели простейшего транспортного средства - двухколесного велосипеда; 2) исследование устойчивости и стабилизация движения неуправляемого двухколесного велосипеда.
Полученные в работе результаты могут быть использованы как базовые при разработке математических моделей, исследовании устойчивости и стабилизации движения четырехколесной техники. Важное значение в этом направлении имеет задача стабилизации движения сельскохозяйственной техники при работе ее на крутых склонах.
Математическая модель движения механизмов типа двухколесного велосипеда представляет также определенный интерес и с точки зрения теоретической механики в плане изучения движения двух или более кинематически связанных колес вокруг мгновенных центров вращения совместно с поступательным движением рамы(платформы).
Данная работа является переработанным и дополненным вариантом работы [1], в которой рассмотрена задача математического моделирования и исследования устойчивости движения двухколесного велосипеда. Первые исследования движения двухколесного велосипеда известны еще с 19 века [2].С момента своего создания велосипеды привлекли внимание многих известных ученых (Список зарубежных работ по данной тематике приведен в конце работы). Публикаций на русском языке значительно меньше. Достаточно полное исследование движения двухколесного велосипеда проведено в монографиях [3,4]. Несмотря на более, чем столетнюю историю плодотворной работы многих ученых, до сих пор не существуетстандартной модели динамики велосипеда и не получены окончательные и исчерпывающие результаты. Тема себя не исчерпала, по ней продолжается работа, появляются новые полезные результаты. Модель движения велосипеда относится к классу задач из области нечеткого моделирования, поэтому имеется достаточно много альтернативных формулировок, различной параметризации (до 25 параметров), различных вариантов динамических переменных, различных способов стабилизации движения. В последнее время интерес к этой проблеме возрос в связи с работами по созданию мобильных роботов [5,6] и возможностью компьютерного моделирования движения велосипеда. В применении к роботам общая постановка задачи математического моделирования может быть сформулирована следующим образом: построить математическую модель опорно-двигательной системы мобильного робота, состоящей из платформы (у велосипеда это рама) и колес, прикрепленных к платформе. Горизонтальные оси колес могут жестко крепиться к платформе (как заднее колесо велосипеда), либо соединяться с платформой с помощью цилиндрического шарнира (как переднее колесо велосипеда). С точки зрения математического моделирования наибольший интерес представляет простейший случай: платформы и двух колес, одно из которых к платформе крепится жестко, а другое шарнирно, что соответствует элементарному велосипеду. Принципиальное отличие методов построения модели в данной работе заключается в применении моделей первого приближения. Модели первого приближения при решении ряда задач позволяют либо получать новые результаты, которые по точным моделям получить не удается, в силу их сложности, либо получать уже известные результаты, но более простым методом. Так в задачах исследования устойчивости движения в случае, когда движение, соответствующее приближенной модели, асимптотически устойчиво, движение, соответствующее точной модели, по теореме Ляпунова об устойчивости по первому приближению, также будет асимптотически устойчивым.
Подобные модели, методы построения моделей, методы исследования устойчивости движения могут использоваться при проектировании различных транспортных средств, их оптимизации, исследования их устойчивого движения, в том числе при расчете условий устойчивого движения сельскохозяйственной техники при движении ее по крутым склонам.
В данной работе рассматривается элементарный велосипед с жесткими колесами, у которого вилка руля имеет вертикальное положение. Это допущение, как показано в [1], не имеет принципиального значения, однако упрощает модель. Полученная в данной работе модель отличается от известных, [3, с. 230], [4, с. 332]. Выведенные на ее основе условия устойчивого движения велосипеда также отличаются от условий, полученных в [3], [4].
Объект исследования. Объектом исследования в данной работе является движение без скольжения по горизонтальной плоскости двухколесного велосипеда с рамой и жесткими колесами одинаковой массы и одинаковыми радиусами.


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



Методы исследования. Основой построения математической модели движения велосипеда являются динамических уравнений Эйлера. Для исследования устойчивости движения, соответствующего построенной математической модели, используется имитационная модель, основанная на критерии Рауса-Гурвица.
Для получения уравнений движения велосипеда вводятся следующие системы координат (рис 2): 1) неподвижная система координат ; 2) подвижная система координат, которая движется поступательно вместе с опорной точкой переднего колеса, (оси системы параллельны соответствующим осям системы ); 3) подвижная система координат, (оси системы параллельны главным осям инерции колеса, однако система не связана с колесом жестко, с системой связана только ось колеса). Аналогично вводятся две системы координат, связанные с задним колесом. , .
При невозмущенном движении система координат совпадает с системой , при этом велосипед движется равномерно и прямолинейно по опорной плоскости . Также движутся и опорные точки и . (Колеса вращаются соответственно вокруг осей ,). Траектория движения велосипеда совпадает с проекцией оси на опорную плоскость. Велосипед движется со скоростью - угловая скорость вращения колеса, - радиус колеса. При возмущенном движении поворот системы координат относительно системы происходит на малые углы в последовательности . Положение заднего колеса в системе определяется соответственно малыми углами . Поскольку указанные углы предполагаются малыми, нет необходимости рассматривать сложные точные математические модели такого движения. Достаточно построить модели первого приближения. Если углы не малые, то преобразование координат при повороте системы координат относительно системы не инвариантно относительно порядка изменения углов. Этот факт не позволяет строить адекватную полную модель движения. В моделях первого приближения преобразование координат инвариантно относительно порядка изменения углов [1]. Мы будем строить модель первого приближения, полагая, что в возмущенном движении угловые скорости поворотов также являются малыми величинами.
В моделях первого приближения будем считать, что для малых величин имеют место соотношения .
Преобразование координат. Последовательные повороты системы координат относительно системы координат на малые углы
, (рис. 2), происходят вместе с колесом и задаются соответственно матрицами [7]

, .
Ориентация системы относительно системы , получаемая в результате двух последовательных поворотов, определяется матрицей

.(1)

Пусть - орты координатных осей , - орты координатных осей . Тогда таблица перехода от векторов к векторам и наоборот с учетом (1) будет иметь вид

1
10
01

Определим далее вектор угловой скорости переднего колеса. В соответствии с рис. 2. в первом приближении будем иметь ,



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

. (2)

Модель движения переднего и заднего колес. Динамические уравнения Эйлера. В работе [5] рекомендуется составлять уравнения движения робота на основе уравнений Лагранжа II - го рода для неголономной системы. Однако, если строить приближенные модели, то предпочтительней использовать динамические уравнения Эйлера, т.к. в этом случае достаточно модели первого приближения. Если использовать уравнения Лагранжа II - го рода, то необходимо строить модель с точностью до малых второго порядка [1].
Не разрывая велосипед на части, (чтобы не возникли реакции связей), составим уравнения движения отдельно для переднего колеса, для заднего и для рамы. Затем эти уравнения объединим.
Для получения уравнений качения переднего колеса воспользуемся теоремой об изменении кинетического момента относительно опорной точки переднего колеса . Пусть - кинетический момент колеса относительно центра (рис. 2), - главный момент внешних сил относительно того же центра. Обозначим массу колеса буквой . Тогда теорема моментов в неподвижной системе координат относительно подвижного в этой системе центра будет иметь вид

, (3)

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

.
.

Выразим кинетический момент в системе .


.

Проекции кинетического момента на оси системы будут

,
.

Тогда

, , (4)
.

Пусть - вес колеса, приложенный в центре масс колеса . Точка имеет в системе координаты Тогда

, ,

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

.

Следовательно,

.(5)

Далее получим выражение . - это скорость перемещения опорной точки вместе с системой при качении колеса по опорной плоскости,

,. ,
.

Следовательно,


.(6)

Поскольку оси системы параллельны соответствующим осям системы , то с учетом (4), (5), (6) в системе координат получаем дифференциальные уравнения движения переднего колеса

. (7)

Это модель качения переднего колеса, построенная на основе динамических уравнений Эйлера в проекциях на неподвижные оси координат.
Кинематическая связь между колесами
Для заднего колеса и для рамы уравнения (3) составим относительно опорной точки заднего колеса .
Модель движения заднего колеса можно получить из модели (7), заменив угол на угол. Угол наклона заднего колеса с точностью до малых первого порядка равен углу наклона переднего колеса [7]. Следовательно, модель движения заднего колеса в системе координат будет иметь вид

. (8)
Получим соотношение между углами и .
3. Кинематическая связь переднего и заднего колес велосипеда.



Рассмотрим траектории движения переднего и заднего колес велосипеда на плоскости . Пусть в некоторый момент времени переднее колесо соприкасается с опорной плоскостью в точке , (рис. 3), заднее в точке . Векторы - соответственно скорости точек соприкосновения переднего и заднего колес, . Звено совершает плоское движение вокруг мгновенного центра , , ,. Угловая скорость звена определяется соотношением

,

где 2а - расстояние между опорными точками переднего и заднего колес.
Обозначим

, тогда


Это уравнение кинематической связи переднего и заднего колес велосипеда.
Рассмотрим случай, когда углы поворота переднего и заднего колес являются малыми. В этом случае . Определим угол поворота переднего колеса формулой , где с,- некоторые константы, [1]. В этом случае траектория движения переднего колеса представляет собой синусоиду.
Тогда последнее уравнение примет вид

.

Решая это линейное дифференциальное уравнение, получим



При имеем.

Тогда .

Таким образом, получим

.

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

, (9)

получим , .
Для достаточно больших можно положить . Тогда .
Уравнения движения рамы. Уравнения движения рамы запишем в неподвижной системе координат относительно подвижного в этой систем центра О. Пусть - кинетический момент рамы относительно центра , - главный момент внешних сил относительно того же центра, - радиус вектор центра масс рамы. Обозначим массу рамы буквой . Тогда теорема моментов в неподвижной системе координат относительно подвижного в этой системе центра будет иметь вид

,

где - количество движения рамы, - скорость центра масс рамы, - угловая скорость вращения рамы вокруг осей и .

.

Кинетический момент рамы выражается соотношением

,

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



Тогда

,
.

В системе имеем ,

,
,
. (10)

Пусть - вес рамы, приложенный в центре масс рамы ,

.

Точка имеет в системе координаты . Тогда

.

где - радиус-вектор центра масс рамы относительно опорной точкиОзаднего колеса. Момент силы относительно центра равен ,

.

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

.

Главный момент внешних сил

(11)


- это скорость перемещения опорной точки вместе с системой при качении колеса по опорной плоскости, ,



. (12)

Поскольку оси системы параллельны соответствующим осям системы , то с учетом (10), (11), (12) в системе координат получаем дифференциальные уравнения движения рамы



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

(13)

Добавим в уравнения движения переднего колеса (7) уравнения движения заднего колеса (8) и уравнения движения рамы (13). Учитывая соотношение , получим полную математическую модель движения велосипеда.

(14)
Исследование устойчивости движения. Исследование устойчивости движения, соответствующего модели (14) с помощью имитационной модели.
Введем обозначения

, ,
, , ,
тогда , , , .

Перепишем уравнения (14) в новых обозначениях

.

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

,
или .

Далее имеем

.

Введем еще обозначения

,


и подставим вместоего выражение (6), получим




Сгруппируем по степеням


.

Далее обозначим


.

Теперь характеристическое уравнение запишется в виде

.

Для исследования устойчивости соответствующего движения применим критерий Рауса-Гурвица [9].
Запишем матрицу из коэффициентов характеристического уравнения в виде


Условие устойчивости:

,



Значение больше нуля. Так как , то при значение также больше нуля. Следовательно, нужно определить значения , при которых и имеют положительные значения.
Исследование значений и проводились численными методами. При этом были приняты следующие значения для величин, входящих в и . м., кг., 2 м., м. Отношение массы рамы к массе колеса задавалось параметром , .
Предполагая, что масса колеса сосредоточена по его ободу, получим моменты инерции колеса . Моменты инерции в формулах (14) определены относительно осей , , . Тогда моменты инерции колеса относительно осей, параллельных осям , , и проходящим через центр масс колеса будут

, , .

Моменты инерции определяются формулами

, ,

. Моменты инерции рамы определялись по упрощенным формулам , (момент инерции стержня относительно оси, проходящей через конец стержня), .
Численные расчеты показали, что больше нуля для любых значений . Это означает, что устойчивость движения велосипеда определяется только знаком определителя
При , когда масса рамы близка к нулю, устойчивое движение велосипеда начинается при угловой скорости рад/с. Если радиус колеса м., то это соответствует линейной скорости велосипеда равной 1 м/с. или 3,6 км/ч. При устойчивое движение велосипеда начинается при угловой скорости рад/с, что соответствует линейной скорости велосипеда равной 2,17 м/с. или 7,8 км/ч. Если, т.е. масса рамы равна массе колеса, минимальная угловая скорость, при которой возможно устойчивое движение велосипеда существенно увеличивается: рад/с, это соответствует линейной скорости велосипеда равной 22,4 м/с. или 80,6 км/ч, (рис. 4). На рис. 5 приведен график зависимости минимальной угловой скорости от параметра при устойчивом движении велосипеда ().
На рис. 6 изображена поверхность, отображающая зависимость от угловой скорости и параметра . При приближении массы рамы к массе колеса устойчивость движения резко снижается.





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

(15)

Обозначения преобразуются к виду

, , .

С учетом этих обозначений уравнения (15) можно записать в виде

(16)

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

,
или .
Далее имеем

Получаем следующий корень характеристического уравнения . Нулевые корни характеристического уравнения свидетельствуют о том, что система (15) допускает частное нулевое решение [7, с. 34].
Подставим вместоего значение (9)

, получим
.
Откуда .

Остальные корни получим из соотношения

, (17)
где , ,

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

.

Следовательно, для устойчивого движения велосипеда угловая скорость должна удовлетворять условию

.(18)
Тогда , где
. (19)

Для велосипеда с приведенными выше параметрами формула (18) дает значение. По критерию Рауса-Гурвица матрица из коэффициентов характеристического уравнения (17) имеет лишь одну ненулевую строку. Все дискриминанты, следовательно, будут нулевыми. При стремлении параметра k к нулю значение стремиться к числу 2,8, что несколько выше значения .
Малые колебания колес велосипеда при его устойчивом движении. Решим систему дифференциальных уравнений (16) и получим выражения, определяющие изменения углов в зависимости от времени . В системе двух уравнений (16) содержится три неизвестные . Значение мы уже определили, следовательно, одной из переменных можно придать произвольное значение, а значение другой вычислить по одной из формул (16).
Пусть. Перенумеруем .

.

Подставив во второе уравнение (16)

, получим
, где

,
,
, где
.

Для из уравнения аналогично получим . Значению будет соответствовать , следовательно,

.

Частные решения системы дифференциальных уравнений (16) будут



,
,


,
, ,
.

За системы частных решений можно отдельно взять действительные и мнимые части.

, ,

Общим решением будет

или
,
, (20)

где константы интегрирования, которые находятся по начальным условиям. Пусть заданы начальные условия: при имеем , .
Определим производные

,
. (21)

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

,, (22)
, (23)
Для угла : .

Тогда углы и будут меняться со временем согласно уравнениям

,
, (24)
где - определены соответственно формулами (19), (22) и (23),
, .

Гироскопическая стабилизация. Как показывает рис. 5 минимальная скорость устойчивого движения велосипеда резко возрастает при стремлении массы рамы к массе колеса. В данном пункте работы исследована устойчивость движения велосипеда при наличии гироскопической стабилизации. Впервые исследования движения твердых тел с гироскопом начали проводиться в конце 19 века. В 1892 г. Бобылев [10] разрешил в эллиптических функциях времени задачу о движении без скольжения по горизонтальной плоскости шара с гироскопом внутри. В то же время эту задачу исследовал Жуковский [11]. Известна классическая задача гироскопической стабилизации однорельсового вагона [3]. В работе [6] рассматривается система стабилизации вертикального положения двухколесного робота-велосипеда.
По теореме Томсона - Тета - Четаева гироскопическую стабилизацию можно осуществить только при четном числе неустойчивых координат. В нашей задаче имеется две неустойчивые координаты: угол θ и угол . Следовательно, гироскопическую стабилизацию движения велосипеда можно осуществить. Пусть в невозмущенном состоянии ось гироскопа, центр масс которого совпадает с центром масс рамы, направлена вдоль оси, проходящей через центр массС рамы параллельно оси . Предположим, что рама повернулась относительно оси на малый угол , а относительно оси на малый угол . Пусть угловая скорость поворота будет . Если гироскоп вращается вокруг собственной оси с угловой скоростью , , то согласно правилу Жуковского со стороны гироскопа будет действовать пара сил с гироскопическим моментом , поворачивающая ось вращения гироскопа в сторону вектора. При этом имеет место соотношение

,

где момент инерции гироскопа.
Добавим в систему уравнений (14) гироскопический момент , получим

(25)

Ведя обозначения

,,
, , ,

где - некоторые константы [1], будем иметь

.

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

,
или .

Подставим вместоего выражение (9), получим



Сгруппируем по степеням и обозначим коэффициенты соответственно при

,,
,

.

Теперь характеристическое уравнение запишется в виде

.

Для исследования устойчивости соответствующего движения применим критерий Рауса-Гурвица. Матрица из коэффициентов характеристического уравнения будет иметь вид



Условие устойчивости:

,



Значение больше нуля. Так как, то при значение также больше нуля. Следовательно, нужно определить значения , при которых и имеют положительные значения.
Исследование значений и проводились численными методами.
Результаты приведены ниже.



Если угловая скорость гироскопа w1=0, то устойчивое движение велосипеда начинается при угловой скорости колесаw>64. При угловой скорости гироскопа w1=1700 область устойчивого движения велосипеда начинается со значения w>1.93



На рис. 8 приведен график зависимости минимальной угловой скорости колеса при устойчивом движении велосипеда с гиростабилизаторомот соотношения массы рамы к массе колеса, угловая скорость вращения гироскопа .
Имитационная модель. Ниже приведена имитационная модель, соответствующая движению велосипеда с гироскопическим стабилизатором. Модель реализована в программной среде MATLAB (программа m21).
Перепишем систему(25), оставив в уравнениях угол поворота заднего колеса и рамы



Введя обозначения

,,
,

перепишем последнюю систему в виде



Выразим из первого уравнения , а из второго уравнения



Обозначим


.

Получим



Массу рамы примем равной массе колеса. Принимая параметры велосипеда из пункта 6, получим значения констант .


Структура модели движения двухколесного велосипеда с гиростабилизатором будет иметь вид



Приведенная модель соответствует значению угловой скорости гироскопа w1=1700. Угловая скорость вращения колеса велосипеда w=2.5. Движение велосипеда устойчивое.


График изменения углав зависимости от времени


Графики изменения углов(желтый) и (красный) в зависимости от времени


График изменения угла в зависимости от времени


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

Выводы

1.Устойчивое движение велосипеда начинается с некоторого критического значения скорости движения велосипеда. При увеличении скорости движения велосипеда устойчивость прямолинейного движения повышается.
.Увеличение соотношения массы рамы к массе колеса ведет к снижению устойчивости движения велосипеда и повышает критическое значение скорости устойчивого движения.
.Устойчивость движения велосипеда можно повысить с помощью гироскопической стабилизации.
Замечание. В данной работе постановка задачи отличается от постановки задачи в работах [3, 4]. В данной работе сразу строится модель первого приближения, а в работах [3, 4] сначала строится точная модель, а затем она линеаризуется. В работах [3, 4] рассматривается статическое положение велосипеда без учета собственного вращения колес. Это приводит к тому, что отличие угла поворота заднего колеса от угла поворота переднего колеса обуславливается наклоном рулевой оси, ее смещением относительно центра масс переднего колеса и величиной так называемого выноса. Это приводит к ошибочному результату. Действительно, если взять элементарный велосипед, у которого рулевая вилка вертикальна, и проходит через центр переднего колеса, то по формулам (6.70), Бутенин [3, с. 235], получается, что угол поворота заднего колеса равен углу поворота переднего колеса. Очевидно, что это не так. Для элементарного велосипеда угол поворота заднего колеса меньше угла поворота переднего колеса [8]. При этом если велосипед не движется, то при любом повороте переднего колеса, заднее колесо остается неподвижным. Если при движении велосипеда угол поворота переднего колеса остается постоянным, то угол поворота заднего колеса изменяется, постепенно приближаясь к углу поворота переднего. Отметим также, что в работе [3] утверждается: неуправляемый велосипед будет двигаться устойчиво в некотором интервале скоростей . В данной работе условие устойчивости движения представлено односторонним неравенством .

Основная литература

1. Бондарь А.И. Построение приближенной модели и исследование устойчивости движения велосипеда. // Труды Новосибирского архитектурно - строительного университета, Т.11, №1, Новосибирск, 2008, с. 78-105, (есть в интернете)
2. Bourlet C. Traite des bicycles etbicyclettes // Paris: Gauthier-Villars, 1898
3. Неймарк Ю.И. Динамика неголономных систем // Ю.И. Неймарк, Н.А. Фуфаев.-М: Наука, 1967. - 519 с.
. Бутенин Н.В., Введение в аналитическую механику. // Н.В. Бутенин, Н.А. Фуфаев. - М., «Наука», 1991, 255 с.
. Мартыненко Ю.Г. Управление движением мобильных роботов., // Мартыненко Ю.Г. Фундаментальная и прикладная математика. - т. 11, вып. 8, 2005, с. 29-80.
. Ленский А.В. Двухколесный робот-велосипед с гиростабилизатором. // А.В. Ленский, А.М. Формальский. Изв. РАН. Теория и сист. управления. - 2003, №3, с. 176-183
. Бондарь А.И. Построение приближенных моделей сферического движения твердого тела // Бондарь И.М., Бондарь А.И - Вычислительные технологии,
Т. 12, №4, 2007, с. 27-41 (есть в интернете)
. Бондарь И.М. Моделирование движения кинематически зависимых тел. // Бондарь И.М., Бондарь А.И - Труды Новосибирского архитектурно - строительного университета, Т.12, №2, Новосибирск, 2009, с. 84-97 (есть в интернете)
. Алфутов Н.А. Устойчивость движения и равновесия // Н.А. Алфутов, К.С. Колесников. - М: Из-во МГТУ им. Баумана, 2003. - 253 с.