Эффективность Спецификаций

Рисунок ниже показывает блок системы с обратной связью. Схема включает номинальную модель самолета, контроллер 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-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.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? The 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 во временном интервале можно сравнить временные характеристики номинальных и наихудших моделей замкнутого цикла с идеальной характеристикой «Handling Quality». Для этого сначала создайте «истинную» модель замкнутой системы 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 рад/сек.