В этом примере показано, как использовать Control System Toolbox™ для разработки цифрового сервопривода для головки чтения/записи дисков.
Для получения дополнительной информации о системе и модели см. Главу 14 «Цифровое управление динамическими системами», Franklin, Powell и Workman.
Ниже приведено изображение системы, которая будет смоделирована.
Узел головного диска (HDA) и приводы моделируются передаточной функцией 10-го порядка, включающей два режима жесткого тела и первые четыре резонанса.
Вход модели является током ic, приводящим в действие двигатель речевой катушки, а выходной сигнал - сигнал ошибки положения (PES, в% от ширины дорожки). Модель также включает небольшую задержку.
Модель дисковода:
Коэффициенты связи, демпфирование и естественные частоты (в Гц) для доминирующих гибких режимов перечислены ниже.
Данные модели:
Учитывая эти данные, создайте номинальную модель сборки головки:
load diskdemo Gr = tf(1e6,[1 12.5 0],'outputdelay',1e-5); Gf1 = tf(w1*[a1 b1*w1],[1 2*z1*w1 w1^2]); % first resonance Gf2 = tf(w2*[a2 b2*w2],[1 2*z2*w2 w2^2]); % second resonance Gf3 = tf(w3*[a3 b3*w3],[1 2*z3*w3 w3^2]); % third resonance Gf4 = tf(w4*[a4 b4*w4],[1 2*z4*w4 w4^2]); % fourth resonance G = Gr * (ss(Gf1) + Gf2 + Gf3 + Gf4); % convert to state space for accuracy
Постройте график отклика Bode модели сборки головки:
cla reset G.InputName = 'ic'; G.OutputName = 'PES'; h = bodeplot(G); title('Bode diagram of the head assembly model'); setoptions(h,'Frequnits','Hz','XLimMode','manual','XLim', {[1 1e5]});
Управление сервоприводом используется, чтобы сохранить головку чтения/записи «на дорожке». Сервопривод контроллера C (z) является цифровым и предназначен для поддержания PES (смещение от центра дорожки) вблизи нуля.
Рассматриваемая здесь нарушение порядка является изменением d шага в вход токе ic. Ваша задача - спроектировать цифровой компенсатор C (z) с адекватной эффективностью подавления помех.
Цифровой сервопривод имеет шаг расчета Ts = 7e-5 с (14,2 кГц).
Спецификации реалистичного проекта перечислены ниже.
Проектные характеристики:
Коэффициент усиления без разомкнутого контура > 20dB при 100 Гц
Пропускная способность > 800 Гц
Запас по амплитуде > 10 дБ
Запас по фазе > 45 o
Пиковое усиление в системе с обратной связью < 4 дБ
Поскольку сервопривод контроллера цифровой, можно выполнить проект в дискретной области. С этого эффекта дискретизируйте модель HDA с помощью C2D и метода удержания нулевого порядка (ZOH):
cla reset Ts = 7e-5; Gd = c2d(G,Ts); h = bodeplot(G,'b',Gd,'r'); % compare with the continuous-time model title('Continuous (blue) and discretized (red) HDA models'); setoptions(h,'Frequnits','Hz','XLimMode','manual','XLim', {[1 1e5]});
Теперь к проекту компенсатора. Начните с чистого интегратора 1/( z-1), чтобы гарантировать нулевую установившуюся ошибку, постройте корневой годограф модели разомкнутого контура Gd * C и масштабируйте z = 1 с помощью опции Zoom In в меню Tools.
C = tf(1,[1 -1],Ts); h = rlocusplot(Gd*C); setoptions(h,'Grid','on','XLimMode','Manual','XLim',{[-1.5,1.5]},... 'YLimMode','Manual','YLim',{[-1,1]});
Из-за двух полюсов в z = 1 цикл сервопривода нестабильен для всех положительных усилений. Чтобы стабилизировать цикл обратной связи, сначала добавьте пару нулей около z = 1.
C = C * zpk([.963,.963],-0.706,1,Ts); h = rlocusplot(Gd*C); setoptions(h,'Grid','on','XLimMode','Manual','XLim',{[-1.25,1.25]},... 'YLimMode','Manual','YLim',{[-1.2,1.2]});
Затем настройте коэффициент усиления цикла, нажав на локус и перетащив черный квадрат внутрь модуля круга. Коэффициент усиления цикла отображается в маркере данных. Коэффициент усиления приблизительно 50 стабилизирует цикл (набор C1 = 50 * C).
C1 = 50 * C;
Теперь симулируйте реакцию с обратной связью на нарушение порядка тока. Нарушение порядка плавно отклоняется, но PES слишком велик (голова отклоняется от центра дорожки на 45% ширины дорожки).
cl_step = feedback(Gd,C1); h = stepplot(cl_step); title('Rejection of a step disturbance (PES = position error)') setoptions(h,'Xlimmode','auto','Ylimmode','auto','Grid','off');
Далее посмотрите на ответ Bode без разомкнутого контура и запасы устойчивости. Коэффициент усиления при 100 Гц составляет всего 15 дБ (по сравнению со спецификацией 20 дБ), и запас по амплитуде только 7dB, поэтому увеличение коэффициента усиления цикла не является опцией.
margin(Gd*C1) diskdemo_aux1(1);
Чтобы освободить место для более высокого низкочастотного усиления, добавьте узкополосный фильтр около резонанса 4000 Гц.
w0 = 4e3 * 2*pi; % notch frequency in rad/sec notch = tf([1 2*0.06*w0 w0^2],[1 2*w0 w0^2]); % continuous-time notch notchd = c2d(notch,Ts,'matched'); % discrete-time notch C2 = C1 * notchd; h = bodeplot(notchd); title('Discrete-time notch filter'); setoptions(h,'FreqUnits','Hz','Grid','on');
Теперь можно безопасно удвоить коэффициент усиления цикла. Результирующие запасы устойчивости и коэффициент усиления при 100 Гц находятся в пределах спецификаций.
C2 = 2 * C2; margin(Gd * C2) diskdemo_aux1(2);
Подавление помех также значительно улучшилось. Теперь PES остается ниже 20% ширины дорожки.
cl_step1 = feedback(Gd,C1); cl_step2 = feedback(Gd,C2); stepplot(cl_step1,'r--',cl_step2,'b') title('2nd-order compensator C1 (red) vs. 4th-order compensator C2 (blue)')
Проверьте, достигнута ли 3dB пиковая характеристика усиления на T = Gd * C/( 1 + Gd * C) (чувствительность к обратной связи):
Gd = c2d(G,Ts); Ts = 7e-5; T = feedback(Gd*C2,1); h = bodeplot(T); title('Peak response of closed-loop sensitivity T(s)') setoptions(h,'PhaseVisible','off','FreqUnits','Hz','Grid','on', ... 'XLimMode','Manual','XLim',{[1e2 1e4]});
Чтобы увидеть пиковое значение, щелкните правой кнопкой мыши по оси и выберите опцию Максимальная чувствительность в меню Characteristics, затем удерживайте мышь над синим маркером или просто щелкните по нему.
Наконец, давайте проанализируем робастность к изменениям в демпфирующих и собственных частотах 2-го и 3-го гибких режимов.
Изменения параметров:
Сгенерируйте массив моделей 16, соответствующих всем комбинациям экстремальных значений z2, w2, z3, w3:
[z2,w2,z3,w3] = ndgrid([.5*z2,1.5*z2],[.9*w2,1.1*w2],[.5*z3,1.5*z3],[.8*w3,1.2*w3]); for j = 1:16, Gf21(:,:,j) = tf(w2(j)*[a2 b2*w2(j)] , [1 2*z2(j)*w2(j) w2(j)^2]); Gf31(:,:,j) = tf(w3(j)*[a3 b3*w3(j)] , [1 2*z3(j)*w3(j) w3(j)^2]); end G1 = Gr * (ss(Gf1) + Gf21 + Gf31 + Gf4);
Дискретизируйте эти 16 моделей сразу и посмотрите, как изменения параметра влияют на характеристику разомкнутого контура. Примечание. Можно кликнуть по любой кривой, чтобы идентифицировать базовую модель.
Gd = c2d(G1,Ts); h = bodeplot(Gd*C2); title('Open-loop response - Monte Carlo analysis') setoptions(h,'XLimMode','manual','XLim',{[8e2 8e3]},'YLimMode','auto',... 'FreqUnits','Hz','MagUnits','dB','PhaseUnits','deg','Grid','on');
Постройте график эффективности подавления помех для этих 16 моделей:
stepplot(feedback(Gd,C2))
title('Step disturbance rejection - Monte Carlo analysis')
Все 16 реакций почти идентичны: наш проект надежен!