Управление оси ответвления самолета Используя Mu-Synthesis

В этом примере показано, как использовать Mu-Analysis и инструменты синтеза в Robust Control Toolbox™. Это описывает проект устойчивого контроллера для боковой направленной оси самолета во время приводимого в действие подхода к приземлению. Линеаризовавшая модель самолета получена для угла нападения 10,5 градусов и скорости полета 140 узлов.

Требования по производительности

Рисунок ниже показов блок-схема системы с обратной связью. Схема включает номинальную модель самолета, контроллер K, а также элементы получая неопределенность модели и цели производительности (см. следующие разделы для деталей).

Рисунок 1: устойчивая система управления для оси ответвления самолета

Цель проекта состоит в том, чтобы заставить самолет эффективно ответить на боковые входные параметры педали палки и руководящего принципа пилота. Требования по производительности включают:

  • Разъединенные ответы от боковой палки p_cmd прокрутить уровень p и от педалей руководящего принципа beta_cmd к углу заноса beta. Боковые педали палки и руководящего принципа имеют максимальное отклонение +/-1 дюйм.

  • Ответ обработки качества (HQ) самолета от боковой палки, чтобы прокрутить уровень p должен совпадать с ответом первого порядка.

HQ_p    = 5.0 * tf(2.0,[1 2.0]);
step(HQ_p), title('Desired response from lateral stick to roll rate (Handling Quality)')

Рисунок 2: Желаемый ответ от боковой палки, чтобы прокрутить уровень.

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

HQ_beta = -2.5 * tf(1.25^2,[1 2.5 1.25^2]);
step(HQ_beta), title('Desired response from rudder pedal to side-slip angle (Handling Quality)')

Рисунок 3: Желаемый ответ от педали руководящего принципа до угла заноса.

  • Приводы стабилизатора имеют +/-20 градусов и +/-пределы на 50 градусов/с на их углу отклонения и уровне отклонения. Приводы руководящего принципа имеют +/-30 градусов и +/-60 угол отклонения градуса/с и ограничения скорости.

  • Три сигнала измерения (прокручивают уровень p, уровень отклонения от курса r, и поперечное ускорение yac ) проникаются фильтры сглаживания второго порядка:

freq = 12.5 * (2*pi);  % 12.5 Hz
zeta = 0.5;
yaw_filt = tf(freq^2,[1 2*zeta*freq freq^2]);
lat_filt = tf(freq^2,[1 2*zeta*freq freq^2]);

freq = 4.1 * (2*pi);  % 4.1 Hz
zeta = 0.7;
roll_filt = tf(freq^2,[1 2*zeta*freq freq^2]);

AAFilters = append(roll_filt,yaw_filt,lat_filt);

От спецификаций до функций взвешивания

Алгоритмы проекта H-бесконечности стремятся минимизировать самое большое усиление с обратной связью через частоту (H-норма-по-бесконечности). Чтобы применить эти инструменты, мы должны сначала переделать технические требования проекта как ограничения на усиления с обратной связью. Мы используем функции взвешивания, чтобы "нормировать" технические требования через частоту и одинаково взвесить каждое требование.

Мы можем выразить спецификации проекта в терминах функций взвешивания можно следующим образом:

  • Чтобы получить пределы на величине отклонения привода и уровне, выберите диагональный, постоянный вес W_act, соответствие стабилизатору и уровню отклонения руководящего принципа и угловым пределам отклонения.

W_act = ss(diag([1/50,1/20,1/60,1/30]));
  • Используйте 3x3 диагональ, фильтр высоких частот W_n смоделировать содержимое частоты шума датчика в уровне списка, уровне отклонения от курса и каналах поперечного ускорения.

W_n = append(0.025,tf(0.0125*[1 1],[1 100]),0.025);
clf, bodemag(W_n(2,2)), title('Sensor noise power as a function of frequency')

Рисунок 4: степень шума Датчика как функция частоты

  • Ответ от ответвления придерживается p и от педали руководящего принципа до beta должен соответствовать качество обработки предназначается для HQ_p и HQ_beta. Это - совпадающая с моделью цель: минимизировать различие (достигают максимума усиление) между желаемыми и фактическими передаточными функциями с обратной связью. Производительность ограничивается из-за нуля правой полуплоскости в модели на уровне 0,002 рад/с, таким образом, точное отслеживание синусоид ниже 0,002 рад/с не возможно. Соответственно, мы взвесим первую качественную спецификацию обработки с полосовым фильтром W_p это подчеркивает частотный диапазон между 0,06 и 30 рад/секунда.

W_p = tf([0.05 2.9 105.93 6.17 0.16],[1 9.19 30.80 18.83 3.95]);
clf, bodemag(W_p), title('Weight on Handling Quality spec')

Рисунок 5: Вес при обработке качественной спецификации.

  • Точно так же выберите W_beta=2*W_p для второй качественной спецификации обработки

W_beta = 2*W_p;

Здесь мы масштабировали веса W_act, W_n\wp, и W_beta так усиление с обратной связью между всеми внешними входными параметрами и все взвесили выходные параметры, меньше 1 на всех частотах.

Номинальная модель самолета

Пилот может управлять боковым направленным ответом самолета с боковыми педалями палки и руководящего принципа. Самолет имеет следующие характеристики:

  • Два входных параметров управления: дифференциальное отклонение стабилизатора delta_stab в градусах, и отклонение руководящего принципа delta_rud в градусах.

  • Три измеренных выходных параметров: уровень списка p в градусе/с, уровень отклонения от курса r в градусе/с и поперечном ускорении yac в g's.

  • Один расчетный выход: угол заноса beta.

Номинальная боковая направленная модель LateralAxis имеет четыре состояния:

  • Боковая скорость v

  • Уровень отклонения от курса r

  • Уровень списка p

  • Угол вращения phi

Эти переменные связаны уравнениями пространства состояний:

x˙=Ax+Bu,y=Cx+Du

где x = [v; r; p; phi], u = [delta_stab; delta_rud], и y = [beta; p; r; yac].

load LateralAxisModel
LateralAxis
LateralAxis =
 
  A = 
               v         r         p       phi
   v      -0.116    -227.3     43.02     31.63
   r     0.00265    -0.259   -0.1445         0
   p    -0.02114    0.6703    -1.365         0
   phi         0    0.1853         1         0
 
  B = 
        delta_stab   delta_rud
   v        0.0622      0.1013
   r     -0.005252    -0.01121
   p      -0.04666    0.003644
   phi           0           0
 
  C = 
                 v          r          p        phi
   beta     0.2469          0          0          0
   p             0          0       57.3          0
   r             0       57.3          0          0
   yac   -0.002827  -0.007877    0.05106          0
 
  D = 
         delta_stab   delta_rud
   beta           0           0
   p              0           0
   r              0           0
   yac     0.002886    0.002273
 
Continuous-time state-space model.

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

A_S = [tf([25 0],[1 25]); tf(25,[1 25])];
A_S.OutputName = {'stab_rate','stab_angle'};

A_R = A_S;
A_R.OutputName = {'rud_rate','rud_angle'};

Составление моделирования ошибок

Номинальная модель только аппроксимирует истинное поведение самолета. Чтобы составлять несмоделированную динамику, можно ввести относительный термин или мультипликативную неопределенность W_in*Delta_G во входе объекта, где ошибочная динамика Delta_G имейте усиление меньше чем 1 через частоты и функцию взвешивания W_in отражает частотные диапазоны, в которых модель более или менее точна. Там обычно больше моделируют ошибки в высоких частотах так W_in высокая передача.

% Normalized error dynamics
Delta_G = ultidyn('Delta_G',[2 2],'Bound',1.0);

% Frequency shaping of error dynamics
w_1 = tf(2.0*[1 4],[1 160]);
w_2 = tf(1.5*[1 20],[1 200]);
W_in = append(w_1,w_2);

bodemag(w_1,'-',w_2,'--')
title('Relative error on nominal model as a function of frequency')
legend('stabilizer','rudder','Location','NorthWest');

Рисунок 6: Относительная погрешность на номинальной модели самолета как функция частоты.

Создавание неопределенной модели динамики самолета

Теперь, когда мы определили количество ошибок моделирования, мы можем создать неопределенную модель динамики самолета, соответствующей пунктирному, окружают рисунок 7 (то же самое как рисунок 1):

Рисунок 7: динамика Самолета.

Используйте connect функционируйте, чтобы объединить номинальную модель LateralAxis корпуса, модели A_S привода и A_R, и описание ошибки моделирования W_in*Delta_G в одну неопределенную модель Plant_unc отображение [delta_stab; delta_rud] на привод и объект выходные параметры:

% Actuator model with modeling uncertainty
Act_unc = append(A_S,A_R) * (eye(2) + W_in*Delta_G);
Act_unc.InputName = {'delta_stab','delta_rud'};

% Nominal aircraft dynamics
Plant_nom = LateralAxis;
Plant_nom.InputName = {'stab_angle','rud_angle'};

% Connect the two subsystems
Inputs = {'delta_stab','delta_rud'};
Outputs = [A_S.y ; A_R.y ; Plant_nom.y];
Plant_unc = connect(Plant_nom,Act_unc,Inputs,Outputs);

Это производит модель Plant_unc неопределенного пространства состояний (USS) из самолета:

Plant_unc
Plant_unc =

  Uncertain continuous-time state-space model with 8 outputs, 2 inputs, 8 states.
  The model uncertainty consists of the following blocks:
    Delta_G: Uncertain 2x2 LTI, peak gain = 1, 1 occurrences

Type "Plant_unc.NominalValue" to see the nominal value, "get(Plant_unc)" to see all properties, and "Plant_unc.Uncertainty" to interact with the uncertain elements.

Анализ, как ошибки моделирования влияют на ответы разомкнутого цикла

Мы можем анализировать эффект моделирования неопределенности путем выбора случайных выборок несмоделированной динамики Delta_G и графический вывод номинальных и встревоженных ответов времени (анализ Монте-Карло). Например, для дифференциального канала стабилизатора, вес неопределенности w_1 подразумевает 5%, моделируя ошибку в низкой частоте, увеличиваясь до 100% после 93 рад/секунда, как подтверждено диаграммой Боде ниже.

% Pick 10 random samples
Plant_unc_sampl = usample(Plant_unc,10);

% Look at response from differential stabilizer to beta
figure('Position',[100,100,560,500]) 
subplot(211), step(Plant_unc.Nominal(5,1),'r+',Plant_unc_sampl(5,1),'b-',10)
legend('Nominal','Perturbed')

subplot(212), bodemag(Plant_unc.Nominal(5,1),'r+',Plant_unc_sampl(5,1),'b-',{0.001,1e3})
legend('Nominal','Perturbed')

Рисунок 8: Переходной процесс и диаграмма Боде.

Разработка контроллера Боковой Оси

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

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

Рисунок 9: модель разомкнутого цикла отображение внешних входных сигналов к связанным с производительностью выходным параметрам.

Чтобы создать эту модель, начните с блок-схемы системы с обратной связью, удалите блок K контроллера и используйте connect вычислить желаемую модель. Как прежде, возможность соединения задана путем маркировки вводов и выводов каждого блока.

Рисунок 10: Блок-схема для того, чтобы создать модель разомкнутого цикла.

% Label block I/Os
AAFilters.u = {'p','r','yac'};    AAFilters.y = 'AAFilt';
W_n.u = 'noise';                  W_n.y = 'Wn';
HQ_p.u = 'p_cmd';                 HQ_p.y = 'HQ_p';
HQ_beta.u = 'beta_cmd';           HQ_beta.y = 'HQ_beta';
W_p.u = 'e_p';                    W_p.y = 'z_p';
W_beta.u = 'e_beta';              W_beta.y = 'z_beta';
W_act.u = [A_S.y ; A_R.y];        W_act.y = 'z_act';

% Specify summing junctions
Sum1 = sumblk('%meas = AAFilt + Wn',{'p_meas','r_meas','yac_meas'});
Sum2 = sumblk('e_p = HQ_p - p');
Sum3 = sumblk('e_beta = HQ_beta - beta');

% Connect everything
OLIC = connect(Plant_unc,AAFilters,W_n,HQ_p,HQ_beta,...
   W_p,W_beta,W_act,Sum1,Sum2,Sum3,...
   {'noise','p_cmd','beta_cmd','delta_stab','delta_rud'},...
   {'z_p','z_beta','z_act','p_cmd','beta_cmd','p_meas','r_meas','yac_meas'});

Это производит неопределенную модель в пространстве состояний

OLIC
OLIC =

  Uncertain continuous-time state-space model with 11 outputs, 7 inputs, 26 states.
  The model uncertainty consists of the following blocks:
    Delta_G: Uncertain 2x2 LTI, peak gain = 1, 1 occurrences

Type "OLIC.NominalValue" to see the nominal value, "get(OLIC)" to see all properties, and "OLIC.Uncertainty" to interact with the uncertain elements.

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

nmeas = 5;		% number of measurements
nctrls = 2;		% number of controls
[kinf,~,gamma_inf] = hinfsyn(OLIC.NominalValue,nmeas,nctrls);
gamma_inf
gamma_inf = 0.9700

Здесь hinfsyn вычисленный контроллер kinf это сохраняет усиление с обратной связью ниже 1, таким образом, спецификациям можно соответствовать для номинальной модели самолета.

Затем выполните Mu-Synthesis, чтобы видеть, можно ли спецификациям соответствовать надежно при принятии во внимание ошибок моделирования (неопределенность Delta_G). Используйте команду musyn выполнять синтез и использовать musynOptions установить сетку частоты, используемую в Mu-Analysis.

fmu = logspace(-2,2,60);
opt = musynOptions('FrequencyGrid',fmu);
[kmu,CLperf] = musyn(OLIC,nmeas,nctrls,opt);
D-K ITERATION SUMMARY:
-----------------------------------------------------------------
                       Robust performance               Fit order
-----------------------------------------------------------------
  Iter         K Step       Peak MU       D Fit             D
    1           5.097        3.487        3.488            12
    2            1.31        1.292        1.312            20
    3           1.243        1.243        1.692            12
    4           1.692        1.543        1.544            16
    5           1.223        1.223         1.55            12
    6           1.534        1.464        1.465            20
    7           1.295        1.293        1.309            12

Best achieved robust performance: 1.22
CLperf
CLperf = 1.2227

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

Сравнение частотного диапазона контроллеров

Сравните производительность и робастность контроллера H-бесконечности kinf и контроллер mu kmu. Вспомните, что спецификации производительности достигаются, когда усиление замкнутого цикла меньше 1 для каждой частоты. Используйте lft функция, чтобы замкнуть круг вокруг каждого контроллера:

clinf = lft(OLIC,kinf);
clmu = lft(OLIC,kmu);

Какова производительность худшего случая (в терминах усиления с обратной связью) каждого контроллера для моделирования ошибок, ограниченных W_in? wcgain команда помогает вам ответить на этот трудный вопрос непосредственно без потребности в обширном gridding и симуляции.

% Compute worst-case gain as a function of frequency
opt = wcOptions('VaryFrequency','on');

% Compute worst-case gain (as a function of frequency) for kinf
[mginf,wcuinf,infoinf] = wcgain(clinf,opt);

% Compute worst-case gain for kmu
[mgmu,wcumu,infomu] = wcgain(clmu,opt);

Можно теперь сравнить номинал и производительность худшего случая для каждого контроллера:

clf
subplot(211)
f = infoinf.Frequency;
gnom = sigma(clinf.NominalValue,f);
semilogx(f,gnom(1,:),'r',f,infoinf.Bounds(:,2),'b');
title('Performance analysis for kinf')
xlabel('Frequency (rad/sec)')
ylabel('Closed-loop gain');
xlim([1e-2 1e2])
legend('Nominal Plant','Worst-Case','Location','NorthWest');

subplot(212)
f = infomu.Frequency;
gnom = sigma(clmu.NominalValue,f);
semilogx(f,gnom(1,:),'r',f,infomu.Bounds(:,2),'b');
title('Performance analysis for kmu')
xlabel('Frequency (rad/sec)')
ylabel('Closed-loop gain');
xlim([1e-2 1e2])
legend('Nominal Plant','Worst-Case','Location','SouthWest');

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

В отличие от этого mu контроллер kmu имеет немного худшую производительность для номинального объекта когда по сравнению с kinf, но это последовательно обеспечивает эту производительность для всех встревоженных моделей (усиление худшего случая около 1.25). mu контроллер поэтому более устойчив к моделированию ошибок.

Валидация временного интервала устойчивого контроллера

Далее протестировать робастность mu контроллера kmu во временном интервале можно сравнить ответы времени номинала и худшего случая модели с обратной связью с идеальной "Обработкой Качества" ответ. Для этого сначала создайте "истинную" модель CLSIM с обратной связью куда все функции взвешивания и образцы модели HQ были удалены:

kmu.u = {'p_cmd','beta_cmd','p_meas','r_meas','yac_meas'};
kmu.y = {'delta_stab','delta_rud'};

AAFilters.y = {'p_meas','r_meas','yac_meas'};

CLSIM = connect(Plant_unc(5:end,:),AAFilters,kmu,{'p_cmd','beta_cmd'},{'p','beta'});

Затем создайте тестовые сигналы u_stick и u_pedal показанный ниже

time = 0:0.02:15;
u_stick = (time>=9 & time<12);
u_pedal = (time>=1 & time<4) - (time>=4 & time<7);

clf
subplot(211), plot(time,u_stick), axis([0 14 -2 2]), title('Lateral stick command')
subplot(212), plot(time,u_pedal), axis([0 14 -2 2]), title('Rudder pedal command')

Можно теперь вычислить и построить идеал, номинал и ответы худшего случая на тестовые команды u_stick и u_pedal.

% Ideal behavior
IdealResp = append(HQ_p,HQ_beta);
IdealResp.y = {'p','beta'};

% Worst-case response
WCResp = usubs(CLSIM,wcumu);

% Compare responses
clf
lsim(IdealResp,'g',CLSIM.NominalValue,'r',WCResp,'b:',[u_stick ; u_pedal],time)
legend('ideal','nominal','perturbed','Location','SouthEast');
title('Closed-loop responses with mu controller KMU')

Ответ с обратной связью почти идентичен для номинала и худшего случая системы с обратной связью. Обратите внимание на то, что ответ уровня списка самолета отслеживает команду уровня списка хорошо первоначально и затем вылетает от этой команды. Это происходит из-за нуля правой полуплоскости в модели самолета в 0,024 рад/секунда.

Смотрите также

|

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте