Как сделать годограф михайлова
Определение устойчивости систем автоматического управления промышленными роботами
Введение
Необходимым условием работоспособности системы автоматического управления (САУ), является её устойчивость. Под устойчивостью принято понимать свойство системы восстанавливать состояние равновесия, из которого она была выведена под влиянием возмущающих факторов после прекращения их воздействия [1].
Постановка задачи
Теория просто и кратко
Анализ устойчивости системы по методу Михайлова сводится к построению характеристического многочлена замкнутой системы (знаменатель передаточной функции), комплексной частотной функции (характеристического вектора):
(1)
где и – соответственно вещественная и мнимая части знаменателя передаточной функции, по виду которой можно судить об устойчивости системы.
Замкнутая САУ устойчива, если комплексная частотная функция , начинаясь на
стрелки начало координат, проходя последовательно n квадрантов, где n – порядок характеристического уравнения системы, т. е.
(2)
Рисунок 1. Амплитудно-фазовые характеристики (годографы) критерия Михайлова: а) – устойчивой системы; б) – неустойчивой системы (1, 2) и системы на границе устойчивости (3)
САУ электроприводом манипулятора промышленного робота (МПР)
Рисунок 2 – Структурная схема САУ электроприводом МПР
Передаточная функция данной САУ имеет следующее выражение [2]:
(3)
где kу – коэффициент усиления усилителя, kм – коэффициент пропорциональности частоты вращения двигателя величине напряжения на якоре, Tу – электромагнитная постоянная времени усилителя, Tм – электромеханическая постоянная времени двигателя с учётом инерции нагрузки (по своим динамическим характеристикам двигатель представляет собой передаточную функцию последовательно соединённых инерционного и интегрирующего звеньев), kдс – коэффициент пропорциональности между входной и выходной величинами датчика скорости, K – коэффициент усиления главной цепи: .
Численные значения в выражение передаточной функции следующие:
K = 100 град / (В∙с); kдс = 0,01 В / (град∙с); Tу = 0,01 с; Tм = 0,1с.
Далее запишем характеристический многочлен замкнутой системы
заменив s на :
(4)
Решение на Python
Здесь следует отметить, что подобные задачи на Python ещё никто не решал, во всяком случае я не нашёл. Это было связано с ограниченными возможностями работы с комплексными числами. С появлением SymPy можно сделать следующее:
Характеристический многочлен замкнутой системы –
-I*T1*T2*w**3 — T1*w**2 — T2*w**2 + I*w + 1
Сразу видим, что многочлен третьей степени. Теперь получим мнимую и действительную части в символьном отображении:
Сразу видим вторую степень действительной части и третью мнимой. Подготовим данные для построения годографа Михайлова. Введём численные значения для T1 и T2, и будем менять частоту от 0 до 100 с шагом 0.1 и построим график:
Из графика не видно, то годограф начинается на действительной положительной оси. Нужно изменить масштабы осей. Приведу полный листинг программы:
Теперь уже видно, что годограф начинается на действительной положительной оси. САУ устойчива, n=3, годограф совпадает с приведённым на первом рисунке.
Дополнительно убедится в том, что годограф начинается на действительной оси можно дополнив программу следующим кодом для w=0:
Начальная точка М(1,0)
САУ сварочного робота
Рисунок 3. Структурная схема САУ позиционированием НСУ
Характеристическое уравнение данной САУ будет иметь вид [1]:
где K – варьируемый коэффициент усиления системы, a – определённая положительная константа. Численные значения: K = 40; a = 0,525.
Далее путём замены s на , получим функцию Михайлова:
(5)
Решение на Python
Построенный годограф Михайлова, начинаясь на вещественной положительной оси (М (21,0)), огибает в положительном направлении начало координат, проходя последовательно четыре квадранта, что соответствует порядку характеристического уравнения. Значит, данная САУ позиционированием НСУ – устойчива.
Выводы
При помощи модуля SymPy Python получен простой и наглядный инструмент для решения задач расчёта устойчивости систем автоматического управления, что является обязательным условием работоспособности любого промышленного робота и манипулятора.
Алгоритм построения годографа Михайлова
Рассмотрим последовательность расчёта критерия устойчивости Михайлова и сформируем алгоритм построения годографа, используя математический пакет «MathCad», на приведенных ниже примерах.
Пример 1. Используя критерий Михайлова, определим устойчивость системы автоматического управления электроприводом манипулятора промышленного робота (МПР). Структурная схема САУ электроприводом МПР изображена на рис. 2.
Рисунок 2 – Структурная схема САУ электроприводом МПР
Передаточная функция данной САУ имеет следующее выражение [2]:
(3)
(4)
где kу – коэффициент усиления усилителя, kм – коэффициент пропорциональности частоты вращения двигателя величине напряжения на якоре, Tу – электромагнитная постоянная времени усилителя, Tм – электромеханическая постоянная времени двигателя с учётом инерции нагрузки (по своим динамическим характеристикам двигатель представляет собой передаточную функцию последовательно соединённых инерционного и интегрирующего звеньев), kдс – коэффициент пропорциональности между входной и выходной величинами датчика скорости, K – коэффициент усиления главной цепи: 
.
Подставим численные значения в выражение передаточной функции:
K = 100 град / (В∙с); kдс = 0,01 В / (град∙с); Tу = 0,01 с; Tм = 0,1с.
(5)
Далее запишем характеристический многочлен замкнутой системы 
заменив s на 
:
(6)
С помощью (1) выделим вещественную и мнимую части и подставим численные значения в полученную комплексную частотную функцию:
(7)
Имея данные в виде (7), перейдём непосредственно к использованию математического пакета «MathCad».
Для этого в верхнем меню выберем «Новый…» – «Пустой документ», в котором будем формировать программу построения годографа Михайлова, используя нижеприведенный алгоритм.
Шаг 1. Задать разрешение годографа диапазоном значений индекса i. Например:
(8)
Шаг 2. Определить исследуемый диапазон и шаг частоты 
, используя значения индекса i (обычно, для практических расчётов, максимальная величина частоты 
не превышает значения 1000, в нашем же примере – достаточно принять 
с частотным шагом 0,1):
(9)
Шаг 3. Полученные вещественную 
и мнимую 
части характеристического уравнения, зададим численными значениями (в данном случае используя (7)) в виде:
(10)
и
(11)
Шаг 4. В результате вычислений (9), (10) и (11), получаются массивы значений частоты 
, а также вещественной 
и мнимой 
частей (рис. 3).
Шаг 5. Далее имея рассчитанные массивы значений 
и 
, переходим к построению годографа Михайлова, используя встроенную функцию «MathCad» – «Инструменты графиков», выбрать «Декартов график». Здесь необходимо определить идентификаторы осей (в данном случае ось абсцисс соответствует вещественной части , а ординат – мнимой части 
) и параметры графика в подменю «Формат…». В результате получим график комплексной частотной функции, приведенный на рис. 4.
Рисунок 4 – Годограф Михайлова для САУ электроприводом МПР
Используя функцию «Трассировка…» (пунктирные линии на рис. 4), можно определить, в соответствующем трассировке окне, точные значения годографа в любой точке рассчитанных массивов.
Таким образом, по рассчитанным данным, построенный годограф Михайлова, начинаясь на действительной положительной полуоси, огибает против часовой стрелки начало координат, проходя последовательно три квадранта, что соответствует порядку характеристического уравнения. Следовательно, данная САУ электроприводом МПР – устойчива.
В соответствии с изложенным алгоритмом, рассмотрим ещё один пример расчёта критерия устойчивости Михайлова и построения комплексной частотной функции.
Пример 2. На современных автомобильных заводах широко применяются большие сварочные роботы (рис. 5). Наконечник сварочного узла (НСУ) подводится к различным местам кузова автомобиля, быстро и точно совершает необходимые действия. Требуется определить устойчивость по критерию Михайлова САУ позиционированием НСУ, структурная схема которой изображена на рис. 6.
Характеристическое уравнение данной САУ будет иметь вид [1]:
(12)
где K – варьируемый коэффициент усиления системы, a – определённая положительная константа.
Подставим в (12) численные значения: K = 40; a = 0,525.
(13)
Далее путём замены s на 
, получим функцию Михайлова:
(14)
Перед тем, как выделить вещественную и мнимую части, запишем (1) в несколько усовершенствованном виде, с целью универсального использования для различных порядков n:
(15)
где с – соответствующий постоянный коэффициент при определённом порядке частоты 
.
Применяя (15) к нашей задаче, получим:
(16)
(17)
Имея данные в виде (16) и (17), приступим к вышеупомянутому алгоритму построения годографа Михайлова с помощью «MathCad».
Шаг 1. Зададим диапазон индекса i:
(18)
Шаг 2. Определим исследуемый диапазон и шаг частоты 
(примем 
с частотным шагом 0,1):
(19)
Шаг 3. Введём вещественную 
(16) и мнимую 
(17) части характеристического уравнения:
(20)
(21)
Шаг 4. При выполнении вычислений (19), (20) и (21), формируются массивы значений частоты 
, вещественной 
и мнимой 
частей (рис. 7).
Шаг 5. Имея рассчитанные массивы значений 
и 
, подобно предыдущему примеру, построим частотную функцию Михайлова. После определения параметров графика, получим годограф, приведенный на рис. 8.
Рисунок 7 – Массивы значений 
, и , рассчитанные в «MathCad»
Рисунок 8 – Годограф Михайлова для САУ позиционированием НСУ
На основе анализа полученных данных, можно сделать вывод, что построенный годограф Михайлова, начинаясь на вещественной положительной оси, огибает в положительном направлении начало координат, проходя последовательно четыре квадранта, что соответствует порядку характеристического уравнения. Значит, данная САУ позиционированием НСУ – устойчива.
Что же касается анализа последних публикаций [2] и сравнения с лучшими аналогами [3], то необходимо отметить факт отсутствия подобной функции (построение годографа Михайлова) в пакете «MATLAB» [1, 4], который, обычно, используется для моделирования различных САУ, что, собственно, и послужило главной причиной создания данного алгоритма.
Перспективы развития данной работы заключаются в создании универсального инструмента для анализа комплексной частотной функции Михайлова, способного выполнить все вычисления уже на этапе задания характеристического уравнения, тем самым полностью автоматизируя этот процесс.
Таким образом, для достижения цели, в ходе написания исследования, была решена главная проблема – получение простого и наглядного инструмента для решения задач расчёта устойчивости САУ, что является обязательным условием работоспособности любого промышленного робота и манипулятора. Также были выполнены следующие задачи: сформирован алгоритм построения комплексной частотной функции Михайлова при помощи математического пакета «MathCad», выполнен анализ устойчивости САУ МПР по данному критерию, кроме того, – приведены практические примеры реализации данного алгоритма.
Список использованной литературы
1. Дорф Р. Современные системы управления / Р. Дорф, Р. Бишоп. – М.: Лаборатория Базовых Знаний, 2002. – 832 с.
2. Юревич Е.И. Основы робототехники 2-е издание / Е.И. Юревич. – С-Пб.: БХВ-Петербург, 2005. – 416 с.
3. Yim Y. Modular Robots / Y. Yim, Y. Zhang, D. Daff // IEEE SPECTRUM. – 2002. – # 2. – P. 30 – 34.
4. Олссон Г. Цифровые системы автоматизации и управления / Г. Олссон, Дж. Пиани. – С-Пб.: Невский Диалект, 2001. – 557 с.
(1)
и 
– соответственно вещественная и мнимая части знаменателя передаточной функции, по виду которой можно судить об устойчивости системы.
и откладывать 
, начинаясь на действительной положительной оси, при изменении частоты 
(2)
(3)
(4)
.
(5)
заменив s на 
:
(6)
(7)
(8)
, используя значения индекса i (обычно, для практических расчётов, максимальная величина частоты 
не превышает значения 1000, в нашем же примере – достаточно принять 
с частотным шагом 0,1):
(9)
части характеристического уравнения, зададим численными значениями (в данном случае используя (7)) в виде:
(10)
и
(11)
, а также вещественной 
и мнимой 
частей (рис. 3).
и 
, переходим к построению годографа Михайлова, используя встроенную функцию «MathCad» – «Инструменты графиков», выбрать «Декартов график». Здесь необходимо определить идентификаторы осей (в данном случае ось абсцисс соответствует вещественной части 
) и параметры графика в подменю «Формат…». В результате получим график комплексной частотной функции, приведенный на рис. 4.
(12)
(13)
, получим функцию Михайлова:
(14)
(15)
.
(16)
(17)
(18)
(примем 
с частотным шагом 0,1):
(19)
(16) и мнимую 
(17) части характеристического уравнения:
(20)
(21)
, вещественной 
и мнимой 
частей (рис. 7).
и 
, подобно предыдущему примеру, построим частотную функцию Михайлова. После определения параметров графика, получим годограф, приведенный на рис. 8.
, 