Устойчивое управление активной приостановки

В этом примере показано, как использовать Robust Control Toolbox™, чтобы спроектировать устойчивый контроллер для активной системы подвески.

Модель приостановки четверти автомобиля

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

Этот пример использует модель четверти автомобиля активной системы подвески (см. рисунок 1). Масса mb (в килограммах), представляет автомобильное шасси (тело) и масса mw (в килограммах), представляет блок колеса. Пружина ks и демпфер bs представляйте пассивную пружину и амортизатор, помещенный между кузовом автомобиля и блоком колеса. Пружина kt моделирует сжимаемость пневматической шины. Переменные xb, xw, и r (все в метрах), перемещение тела, перемещение колеса и дорожное воздействие, соответственно. Сила fs (в килоньютонах), примененный между телом и блоком колеса управляется обратной связью и представляет активный компонент системы подвески.

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

С обозначением (x1,x2,x3,x4):=(xb,xb˙,xw,x˙w), линеаризовавшие уравнения пространства состояний для модели четверти автомобиля:

x1˙=x2x2˙=-(1/mb)[ks(x1-x3)+bs(x2-x4)-103fs]x3˙=x4x4˙=(1/mw)[ks(x1-x3)+bs(x2-x4)-kt(x3-r)-103fs].

Создайте модель в пространстве состояний qcar представление этих уравнений.

% Physical parameters
mb = 300;    % kg
mw = 60;     % kg
bs = 1000;   % N/m/s
ks = 16000 ; % N/m
kt = 190000; % N/m

% State matrices
A = [ 0 1 0 0; [-ks -bs ks bs]/mb ; ...
      0 0 0 1; [ks bs -ks-kt -bs]/mw];
B = [ 0 0; 0 1e3/mb ; 0 0 ; [kt -1e3]/mw];
C = [1 0 0 0; 1 0 -1 0; A(2,:)];
D = [0 0; 0 0; B(2,:)];

qcar = ss(A,B,C,D);
qcar.StateName = {'body travel (m)';'body vel (m/s)';...
          'wheel travel (m)';'wheel vel (m/s)'};
qcar.InputName = {'r';'fs'};
qcar.OutputName = {'xb';'sd';'ab'};

Передаточная функция от привода до перемещения тела и ускорения имеет нуль мнимой оси с собственной частотой 56,27 рад/с. Это называется частотой транзитного участка шины.

tzero(qcar({'xb','ab'},'fs'))
ans = 2×1 complex

  -0.0000 +56.2731i
  -0.0000 -56.2731i

Точно так же передаточная функция от привода до отклонения приостановки имеет нуль мнимой оси с собственной частотой 22,97 рад/с. Это называется rattlespace частотой.

zero(qcar('sd','fs'))
ans = 2×1 complex

   0.0000 +22.9734i
   0.0000 -22.9734i

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

bodemag(qcar({'ab','sd'},'r'),'b',qcar({'ab','sd'},'fs'),'r',{1 100});
legend('Road disturbance (r)','Actuator force (fs)','location','SouthWest')
title({'Gain from road dist (r) and actuator force (fs) ';
       'to body accel (ab) and suspension travel (sd)'})

Из-за нулей мнимой оси управление с обратной связью не может улучшить ответ от дорожного воздействия r к ускорению тела ab на частоте транзитного участка шины, и от r к отклонению приостановки sd на rattlespace частоте. Кроме того, из-за отношения xw=xb-sd и то, что положение колеса xw примерно следует r в низкой частоте (меньше чем 5 рад/с) существует свойственный компромисс между пассажирским комфортом и отклонением приостановки: любое сокращение перемещения тела в низкой частоте приведет к увеличению отклонения приостановки.

Неопределенная модель привода

Гидравлический привод, используемый в активном управлении приостановкой, соединяется между массой тела mb и масса блока колеса mw. Номинальные движущие силы привода представлены передаточной функцией первого порядка 1/(1+s/60) смещение имеющее 0,05 м.

ActNom = tf(1,[1/60 1]);

Эта номинальная модель только аппроксимирует физическую динамику привода. Мы можем использовать семейство моделей привода, чтобы составлять моделирование ошибок и изменчивости в моделях четвертей автомобиля и приводе. Это семейство состоит из номинальной модели с зависимой частотой суммой неопределенности. В низкой частоте, ниже 3 рад/с, модель может варьироваться до 40% от своей номинальной стоимости. Приблизительно 3 рад/с, изменение процента начинает увеличиваться. Неопределенность пересекает 100% на уровне 15 рад/с и достигает 2 000% на уровне приблизительно 1 000 рад/с. Функция взвешивания Wunc используется, чтобы модулировать сумму неопределенности с частотой.

Wunc = makeweight(0.40,15,3);
unc = ultidyn('unc',[1 1],'SampleStateDim',5);
Act = ActNom*(1 + Wunc*unc);
Act.InputName = 'u';
Act.OutputName = 'fs';

Результат Act неопределенная модель в пространстве состояний привода. Постройте Предвещать ответ 20 демонстрационных значений Act и сравните с номинальной стоимостью.

rng('default')
bode(Act,'b',Act.NominalValue,'r+',logspace(-1,3,120))

Спроектируйте Setup

Основные цели управления формулируются в терминах пассажирского комфорта и дорожной обработки, которые относятся к ускорению тела ab и перемещение приостановки sd. Другие факторы, которые влияют на систему управления, включают характеристики дорожного воздействия, качество измерений датчика для обратной связи и пределы на силе доступного элемента управления. Использовать H алгоритмы синтеза, мы должны выразить эти цели как одну функцию стоимости, которая будет минимизирована. Это может быть сделано как обозначенный рисунок 2.

Рисунок 2: формулировка Подавления помех.

Диспетчер обратной связи использует измерения y1,y2 из перемещения приостановки sd и ускорение тела ab вычислить управляющий сигнал u управление гидравлическим приводом. Существует три внешних источника воздействия:

  • Дорожное воздействие r, смоделированный как нормированный сигнал d1 сформированный функцией взвешивания Wroad. К широкополосным дорожным отклонениям модели величины семь сантиметров мы используем постоянный вес Wroad=0.07

  • Шум датчика на обоих измерениях, смоделированных как нормированные сигналы d2 и d3 сформированный функциями взвешивания Wd2 и Wd3. Мы используем Wd2=0.01 и Wd3=0.5 к широкополосному шуму датчика модели интенсивности 0.01 и 0.5, соответственно. В более реалистическом проекте эти веса были бы зависимым частоты, чтобы смоделировать шумовой спектр ускоряющих датчиков и смещения.

Целям управления можно дать иное толкование как цель подавления помех: Минимизируйте удар воздействий d1,d2,d3 на взвешенной комбинации усилия по управлению u, перемещение приостановки sd, и ускорение тела ab. При использовании H норма (пиковое усиление), чтобы измерить "удар", это составляет разработку контроллера, который минимизирует H норма от входных параметров воздействия d1,d2,d3 к сигналам ошибки e1,e2,e3.

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

Wroad = ss(0.07);  Wroad.u = 'd1';   Wroad.y = 'r';
Wact = 0.8*tf([1 50],[1 500]);  Wact.u = 'u';  Wact.y = 'e1';
Wd2 = ss(0.01);  Wd2.u = 'd2';   Wd2.y = 'Wd2';
Wd3 = ss(0.5);   Wd3.u = 'd3';   Wd3.y = 'Wd3';

Задайте цели с обратной связью для усиления от дорожного воздействия r к отклонению приостановки sd (обработка) и ускорение тела ab (комфорт). Из-за неопределенности привода и нулей мнимой оси, только стремитесь ослабить воздействия ниже 10 рад/с.

HandlingTarget = 0.04 * tf([1/8 1],[1/80 1]);
ComfortTarget = 0.4 * tf([1/0.45 1],[1/150 1]);

Targets = [HandlingTarget ; ComfortTarget];
bodemag(qcar({'sd','ab'},'r')*Wroad,'b',Targets,'r--',{1,1000}), grid
title('Response to road disturbance')
legend('Open-loop','Closed-loop target')

Соответствующие веса производительности Wsd,Wab обратные величины их, успокаивают и обрабатывающие цели. Чтобы исследовать компромисс между пассажирским комфортом и дорожной обработкой, создайте три набора весов (βWsd,(1-β)Wab) соответствие трем различным компромиссам: комфорт (β=0.01), сбалансированный (β=0.5), и обработка (β=0.99).

% Three design points
beta = reshape([0.01 0.5 0.99],[1 1 3]);
Wsd = beta / HandlingTarget;
Wsd.u = 'sd';  Wsd.y = 'e3';
Wab = (1-beta) / ComfortTarget;
Wab.u = 'ab';  Wab.y = 'e2';

Наконец, используйте connect создать модель qcaric из блок-схемы рисунка 2. Обратите внимание на то, что qcaric массив трех моделей, один для каждой точки проекта β. Кроме того, qcaric неопределенная модель, поскольку она содержит неопределенную модель Act привода.

sdmeas  = sumblk('y1 = sd+Wd2');
abmeas = sumblk('y2 = ab+Wd3');
ICinputs = {'d1';'d2';'d3';'u'};
ICoutputs = {'e1';'e2';'e3';'y1';'y2'};
qcaric = connect(qcar(2:3,:),Act,Wroad,Wact,Wab,Wsd,Wd2,Wd3,...
                 sdmeas,abmeas,ICinputs,ICoutputs)
qcaric =

  3x1 array of uncertain continuous-time state-space models.
  Each model has 5 outputs, 4 inputs, 9 states, and the following uncertain blocks:
    unc: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences

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

Номинальный проект H-бесконечности

Используйте hinfsyn вычислить H контроллер для каждого значения смешивающегося фактора β.

ncont = 1; % one control signal, u
nmeas = 2; % two measurement signals, sd and ab
K = ss(zeros(ncont,nmeas,3));
gamma = zeros(3,1);
for i=1:3
   [K(:,:,i),~,gamma(i)] = hinfsyn(qcaric(:,:,i),nmeas,ncont);
end

gamma
gamma = 3×1

    0.9405
    0.6727
    0.8892

Эти три контроллера достигают с обратной связью H нормы 0,94, 0.67 и 0.89, соответственно. Создайте соответствующие модели с обратной связью и сравните усиления от дорожного воздействия до xb,sd,ab для пассивных и активных приостановок. Заметьте, что все три контроллера уменьшают отклонение приостановки и ускорение тела ниже rattlespace частоты (23 рад/с).

% Closed-loop models
K.u = {'sd','ab'};  K.y = 'u';
CL = connect(qcar,Act.Nominal,K,'r',{'xb';'sd';'ab'});

bodemag(qcar(:,'r'),'b', CL(:,:,1),'r-.', ...
   CL(:,:,2),'m-.', CL(:,:,3),'k-.',{1,140}), grid
legend('Open-loop','Comfort','Balanced','Handling','location','SouthEast')
title('Body travel, suspension deflection, and body acceleration due to road')

Оценка временного интервала

Чтобы далее оценить три проекта, выполните симуляции временного интервала с помощью дорожного сигнала воздействия r(t) представление дорожного удара высоты 5 см.

% Road disturbance
t = 0:0.0025:1;
roaddist = zeros(size(t));
roaddist(1:101) = 0.025*(1-cos(8*pi*t(1:101)));

% Closed-loop model
SIMK = connect(qcar,Act.Nominal,K,'r',{'xb';'sd';'ab';'fs'});

% Simulate
p1 = lsim(qcar(:,1),roaddist,t);
y1 = lsim(SIMK(1:4,1,1),roaddist,t);
y2 = lsim(SIMK(1:4,1,2),roaddist,t);
y3 = lsim(SIMK(1:4,1,3),roaddist,t);

% Plot results
subplot(211)
plot(t,p1(:,1),'b',t,y1(:,1),'r.',t,y2(:,1),'m.',t,y3(:,1),'k.',t,roaddist,'g')
title('Body travel'), ylabel('x_b (m)')
subplot(212)
plot(t,p1(:,3),'b',t,y1(:,3),'r.',t,y2(:,3),'m.',t,y3(:,3),'k.',t,roaddist,'g')
title('Body acceleration'), ylabel('a_b (m/s^2)')

subplot(211)
plot(t,p1(:,2),'b',t,y1(:,2),'r.',t,y2(:,2),'m.',t,y3(:,2),'k.',t,roaddist,'g')
title('Suspension deflection'), xlabel('Time (s)'), ylabel('s_d (m)')
subplot(212)
plot(t,zeros(size(t)),'b',t,y1(:,4),'r.',t,y2(:,4),'m.',t,y3(:,4),'k.',t,roaddist,'g')
title('Control force'), xlabel('Time (s)'), ylabel('f_s (kN)')
legend('Open-loop','Comfort','Balanced','Handling','Road Disturbance','location','SouthEast')

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

Устойчивый проект Му

До сих пор вы спроектировали H контроллеры, которые достигают целей производительности для номинальной модели привода. Как обсуждено ранее, эта модель является только приближением истинного привода, и необходимо убедиться, что производительность контроллера обеспечена перед лицом ошибок модели и неопределенности. Это называется устойчивой производительностью.

Затем используйте μ- синтез, чтобы спроектировать контроллер, который достигает устойчивой производительности для целого семейства моделей привода. Устойчивый контроллер синтезируется с функцией musyn использование неопределенной модели qcaric(:,:,2) соответствие "сбалансированной" производительности (β=0.5).

[Krob,rpMU] = musyn(qcaric(:,:,2),nmeas,ncont);
D-K ITERATION SUMMARY:
-----------------------------------------------------------------
                    Closed-loop performance             Fit order
-----------------------------------------------------------------
  Iter         K-Step       Peak MU      D-Step             D
    1           1.193        1.125        1.139             4
    2           1.091        1.025        1.033             4
    3          0.9991        0.946       0.9559             4
    4          0.9358       0.9319       0.9348             4
    5          0.9096       0.9057       0.9114             8
    6          0.9103        0.907       0.9096             8
    7          0.9091       0.9066       0.9094             6

Best achieved robust performance: 0.906

Симулируйте номинальный ответ на дорожный удар с устойчивым контроллером Krob. Ответы похожи на полученных со "сбалансированным" H контроллер.

% Closed-loop model (nominal)
Krob.u = {'sd','ab'};
Krob.y = 'u';
SIMKrob = connect(qcar,Act.Nominal,Krob,'r',{'xb';'sd';'ab';'fs'});

% Simulate
p1 = lsim(qcar(:,1),roaddist,t);
y1 = lsim(SIMKrob(1:4,1),roaddist,t);

% Plot results
clf, subplot(221)
plot(t,p1(:,1),'b',t,y1(:,1),'r',t,roaddist,'g')
title('Body travel'), ylabel('x_b (m)')
subplot(222)
plot(t,p1(:,3),'b',t,y1(:,3),'r')
title('Body acceleration'), ylabel('a_b (m/s^2)')
subplot(223)
plot(t,p1(:,2),'b',t,y1(:,2),'r')
title('Suspension deflection'), xlabel('Time (s)'), ylabel('s_d (m)')
subplot(224)
plot(t,zeros(size(t)),'b',t,y1(:,4),'r')
title('Control force'), xlabel('Time (s)'), ylabel('f_s (kN)')
legend('Open-loop','Robust design','location','SouthEast')

Затем симулируйте ответ на дорожный удар для 100 моделей привода, случайным образом выбранных из неопределенного набора модели Act.

rng('default'), nsamp = 100;  clf

% Uncertain closed-loop model with balanced H-infinity controller
CLU = connect(qcar,Act,K(:,:,2),'r',{'xb','sd','ab'});
lsim(usample(CLU,nsamp),'b',CLU.Nominal,'r',roaddist,t)
title('Nominal "balanced" design')
legend('Perturbed','Nominal','location','SouthEast')

% Uncertain closed-loop model with balanced robust controller
CLU = connect(qcar,Act,Krob,'r',{'xb','sd','ab'});
lsim(usample(CLU,nsamp),'b',CLU.Nominal,'r',roaddist,t)
title('Robust "balanced" design')
legend('Perturbed','Nominal','location','SouthEast')

Устойчивый контроллер Krob уменьшает изменчивость из-за неопределенности модели и поставляет более сопоставимую производительность.

Упрощение контроллера: сокращение порядка

Устойчивый контроллер Krob имеет относительно старший разряд по сравнению с объектом. Можно использовать функции снижения сложности модели, чтобы найти контроллер более низкоуровневый, который достигает того же уровня устойчивой производительности. Используйте reduce сгенерировать приближения различных порядков.

% Create array of reduced-order controllers
NS = order(Krob);
StateOrders = 1:NS;
Kred = reduce(Krob,StateOrders);

Затем используйте robgain вычислить устойчивое поле производительности для каждого приближения уменьшаемого порядка. Целям производительности удовлетворяют, когда усиление с обратной связью меньше γ=1. Устойчивое поле производительности измеряется, сколько неопределенности может быть поддержано без ухудшающейся производительности (превышение γ=1). Поле 1 или больше указывает, что мы можем выдержать 100% заданной неопределенности.

% Compute robust performance margin for each reduced controller
gamma = 1;
CLP = lft(qcaric(:,:,2),Kred);
for k=1:NS
   PM(k) = robgain(CLP(:,:,k),gamma);
end

% Compare robust performance of reduced- and full-order controllers
PMfull = PM(end).LowerBound;
plot(StateOrders,[PM.LowerBound],'b-o',...
   StateOrders,repmat(PMfull,[1 NS]),'r');
grid
title('Robust performance margin as a function of controller order')
legend('Reduced order','Full order','location','SouthEast')

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

Упрощение контроллера: настройка Фиксированного Порядка

В качестве альтернативы можно использовать musyn непосредственно настроить контроллеры младшего разряда. Это часто более эффективно, чем по опыту сокращение контроллера полного порядка Krob. Например, настройте контроллер третьего порядка, чтобы оптимизировать его устойчивую производительность.

% Create tunable 3rd-order controller 
K = tunableSS('K',3,ncont,nmeas);

% Tune robust performance of closed-loop system CL
CL0 = lft(qcaric(:,:,2),K);
[CL,RP] = musyn(CL0);
D-K ITERATION SUMMARY:
-----------------------------------------------------------------
                    Closed-loop performance             Fit order
-----------------------------------------------------------------
  Iter         K-Step       Peak MU      D-Step             D
    1           1.189        1.104        1.121            10
    2           1.077       0.9998        1.014            10
    3          0.9885        0.983       0.9941            10
    4          0.9308       0.9294       0.9378            10
    5          0.9258        0.924       0.9271            10
    6          0.9256       0.9226       0.9364            10

Best achieved robust performance: 0.923

У настроенного контроллера есть производительность RP=0.92, очень близко к тому из Krob. Вы видите, что Предвещать использование ответа

K3 = getBlockValue(CL,'K');
bode(K3)