Технические характеристики

На рисунке ниже показана блок-схема системы с замкнутым контуром. Схема включает номинальную модель самолета, контроллер 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);

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

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-синтез, чтобы увидеть, могут ли спецификации быть выполнены надежно при учете ошибок моделирования (неопределенность Delta_G). Используйте команду musyn для выполнения синтеза и использования musynOptions для установки частотной сетки, используемой для анализа mu.

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.242        1.242        1.693            12
    4           1.693        1.544        1.545            16
    5           1.223        1.223        1.551            12
    6           1.533        1.464        1.465            20
    7           1.278        1.277        1.302            12

Best achieved robust performance: 1.22

CLperf

CLperf = 1.2226

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

Сравнение контроллеров в частотной области

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

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

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

% 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 рад/сек.