Моделирование линейных цифровых систем автоматического управления в MATLAB

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

(18)

где J - момент инерции рычага, на котором установлен блок магнитной головки;

С - коэффициент вязкого трения в опоре рычага;

K – коэффициент жесткости пружины;

K_i – коэффициент передачи двигателя (по моменту);

i (t) – ток якоря двигателя;

j (t) – угловое положение рычага головки.

Если принять за входную величину ток двигателя i (t), а за выходную – угол поворота j (t), то, применяя к уравнению (1) преобразование Лапласа, получим

. (19)

Создадим tf -модель при следующих исходных данных: J = 0,01 кг×м², K = 10 Н×м/рад, C = 0,004 Н×м/(рад/с), K_i = 0,05 Н×м/А:

» J = 0.01; C = 0.004; K = 10; Ki = 0.05;

» num = Ki;

» den = [J C K];

» H = tf(num,den)

Transfer function:

0.05

-----------------------

0.01s^2 + 0.004s + 10

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

» Ts = 0.005;

» Hd = c2d(H,Ts,'zoh')

Transfer function:

6.233e-005z + 6.229e-005

-------------------------

z^2 - 1.973z + 0.998

Sampling time: 0.005

Для преобразования непрерывной модели в дискретную здесь использована функция c2d, у которой период квантования Ts = 0.005 c, а для дискретизации выбран метод 'zoh'. Проверку качества выполненной дискретизации непрерывного объекта можно выполнить, рассмотрев частотные характеристики непрерывной и дискретной моделей объектов (рис. 20):

» bode(H,'-',Hd,'--')

Если в дискретной системе автоматического позиционирования магнитной головки жесткого диска не использовать корректирующие устройства, то переходный процесс, определяемый функцией step, имеет вид(рис. 21):

» step(Hd)

Рис. 20. Частотные характеристики объекта	Рис. 21. Переходная характеристика объекта

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

» damp(Hd);

Eigenvalue	Magnitude	Equiv.Damping	Equiv. Freq.(rad/s)
9.87e-001+1.57e-001i	9.99e-001	6.32e-003	3.16e+001
9.87e-001-1.57e-001i	9.99e-001	6.32e-003	3.16e+001

Отсюда видно, что полюса расположены вблизи окружности единичного радиуса и система, вследствие этого, имеет слабое демпфирование. Таким образом, необходим дискретный регулятор, который увеличит демпфирование в системе. Решим поставленную задачу с помощью метода корневого годографа (МКГ) (см. лабораторную работу № 1). Для этого с помощью функции rltool выполним загрузку основного окна интерфейса МКГ. Далее из меню File с помощью Import Model импортируем модель Hd в блок объекта G. В результате получим распределение (рис. 22) нулей и полюсов модели.

Рис. 22. Распределение нулей и полюсов модели объекта	Рис. 23. Распределение нулей и полюсов модели системы управления

Если выбрать параметры компенсатора с передаточной функцией, показанной на рис. 23, то получим переходный процесс в системе, изображенный на рис. 24.

Рис. 24. Переходный процесс в системе с компенсатором	Рис. 25. Частотные характеристики исходной и скорректированной систем

Если параметры полученного переходного процесса не отвечают требованиям на проектирование, то, изменяя положение полюса и значение коэффициента передачи компенсатора, получают оптимальный переходный процесс. В частности, полученное быстродействие системы (около 0,1 с) не отвечает современным требованиям к управлению магнитными головками. Можно попытаться усложнить передаточную функцию компенсатора, добавив нули и полюса, хотя демпфирование магнитной головки низкое и добиться высокого быстродействия затруднительно.

Синтез дискретного корректирующего устройства рассматриваемой системы можно выполнить с помощью функций Control System Toolbox. В частности после определения дискретной передаточной функции объекта дискретная модель компенсатора может быть получена с помощью функции

» Komp=tf([1 -0.85],[1 0],Ts)%Формирование модели компенсатора

Transfer function:

z - 0.85

--------

z

Sampling time: 0.005

» Hp=Hd*Komp %Формирование разомкнутого контура

С помощью функции

» bode(Hd,'--',Hp,'-')

можно выяснить, как модифицирует выбранный компенсатор частотные характеристики системы (рис. 25).

	Анализируя частотные характеристики системы, видим, что компенсатор сдвигает фазовую характеристику вверх (вводится опережение по фазе). Корневой годограф контура с компенсатором (рис. 26) можно определить с помощью функции » rlocus(Hp), zgrid Далее используя функцию » k,poles]=rlocfind(Hp), задают оптимальное расположение полюсов:
Рис. 26. Корневой годограф контура с компенсатором

Select a point in the graphics window

selected_point =

0.2995 - 0.1024i

k = 3.5241e+003

poles =

0.7081 + 0.2276i

0.7081 - 0.2276i

0.3373

Определим демпфирование замкнутой системы:

» ddamp(poles,Ts)

Eigenvalue	Magnitude	Equiv.Damping	Equiv. Freq.(rad/s)
7.08e-001+2.28e-001i	7.44e-001	6.90e-001	8.59e+001
7.08e-001-2.28e-001i	7.44e-001	6.90e-001	8.59e+001
3.37e-001	3.37e-001	1.00e+000	2.17e+002

Чтобы оценить результаты синтеза системы, сформируем замкнутую систему и построим для нее переходную функцию:

» Hsam=feedback(Hp,k)

Transfer function:

6.233e-005 z^2 + 9.308e-006 z - 5.294e-005

------------------------------------------

z^3 - 1.753 z^2 + 1.031 z - 0.1866

Sampling time: 0.005

» step(Hsam)

Для определения показателей устойчивости системы воспользуемся функцией margin, предварительно учтя в передаточной функции разомкнутой системы найденной с помощью функции rlocfind значение k:

» olk=k*Hp;

» [Gm,Pm,Wcg,Wcp]=margin(olk);

» Margins=[Gm Wcg Pm Wcp]

Margins =

4.4302 296.7989 44.5498 94.8787