Документация

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

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

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

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

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: Желаемый ответ от боковой палки, чтобы прокрутить уровень.

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: Желаемый ответ от педали руководящего принципа до угла заноса.

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 = ss(diag([1/50,1/20,1/60,1/30]));

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: степень шума Датчика как функция частоты

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_act, W_n, W_p и W_beta так усиление с обратной связью между всеми внешними входными параметрами, и все взвесили выходные параметры, меньше чем 1 на всех частотах.

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

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

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

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

где 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). Используйте команду dksyn, чтобы выполнить синтез и использовать dksynOptions, чтобы установить сетку частоты, используемую для Mu-Analysis и количества итераций D-K.

fmu = logspace(-2,2,60);
opt = dksynOptions('FrequencyVector',fmu,'NumberofAutoIterations',5);
[kmu,~,gamma_mu] = dksyn(OLIC,nmeas,nctrls,opt);
gamma_mu

gamma_mu = 1.2228

Здесь лучший 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 радах/секунда.