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

В этом первом примере у нас есть семейство моделей, описывающих поведение объекта в различных рабочих условиях. Номинальная модель объекта управления является нестабильной системой первого порядка.

Pnom = tf(2,[1 -2])

Pnom =
 
    2
  -----
  s - 2
 
Continuous-time transfer function.

Другие модели являются изменениями Pnom. Все они имеют один нестабильный полюс, но расположение этого полюса может варьироваться в зависимости от условий работы.

p1 = Pnom*tf(1,[.06 1]);              % extra lag
p2 = Pnom*tf([-.02 1],[.02 1]);       % time delay
p3 = Pnom*tf(50^2,[1 2*.1*50 50^2]);  % high frequency resonance
p4 = Pnom*tf(70^2,[1 2*.2*70 70^2]);  % high frequency resonance
p5 = tf(2.4,[1 -2.2]);                % pole/gain migration
p6 = tf(1.6,[1 -1.8]);                % pole/gain migration

Чтобы применить устойчивые инструменты управления, мы можем заменить этот набор моделей одной неопределенной моделью объекта управления, область значений поведения включает p1 через p6. Это одно из применений команды ucover. Эта команда принимает массив моделей LTI Parray и номинальную модель Pnom и моделирует различие Parray-Pnom как мультипликативная неопределенность в динамике системы.

Потому что ucover ожидает массив моделей, используйте stack команда для сбора моделей объекта управления p1 через p6 в один массив.

Parray = stack(1,p1,p2,p3,p4,p5,p6);

Далее используйте ucover чтобы «охватить» область значений поведения Parray с неопределенной моделью вида

P = Pnom * (1 + Wt * Delta)

где вся неопределенность сосредоточена в «неизмененной динамике» Delta (а ultidyn объект). Потому что выигрыш Delta равномерно ограничен 1 на всех частотах, «формирующий» фильтр Wt используется для захвата того, как относительное количество неопределенности изменяется с частотой. Этот фильтр также упоминается как функция взвешивания неопределенности.

Попробуйте фильтр 4-го порядка Wt для этого примера:

orderWt = 4;
Parrayg = frd(Parray,logspace(-1,3,60));
[P,Info] = ucover(Parrayg,Pnom,orderWt,'InputMult');

Получившаяся модель P является объектом с одним входом, одним выходом с неопределенным пространством состояний (USS) с номинальным значением Pnom.

P

P =

  Uncertain continuous-time state-space model with 1 outputs, 1 inputs, 5 states.
  The model uncertainty consists of the following blocks:
    Parrayg_InputMultDelta: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences

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

tf(P.NominalValue)

ans =
 
    2
  -----
  s - 2
 
Continuous-time transfer function.

График величины Бода подтверждает, что формирующий фильтр Wt «покрывает» относительные изменения в поведении объектов. Как функция от частоты, уровень неопределенности составляет 30% при 5 рад/сек (-10dB = 0,3), 50% при 10 рад/сек и 100% после 29 рад/сек.

Wt = Info.W1;
bodemag((Pnom-Parray)/Pnom,'b--',Wt,'r'); grid
title('Relative Gaps vs. Magnitude of Wt')

Теперь можно использовать неопределенную модель P для разработки робастного контроллера для исходного семейства моделей объекта управления, см. «Одновременная стабилизация с использованием робастного управления» для получения дополнительной информации.

Упрощение существующей неопределенной модели

В этом втором примере мы начинаем с детальной неопределенной модели объекта. Эта модель состоит из динамики первого порядка с неопределенным коэффициентом усиления и постоянным временем, последовательно с умеренно заниженным резонансом и значительной немодулированной динамикой. Эта модель создается с помощью ureal и ultidyn команды для определения неопределенных переменных:

gamma = ureal('gamma',2,'Perc',30);  % uncertain gain
tau = ureal('tau',1,'Perc',30);      % uncertain time-constant
wn = 50; xi = 0.25;
P = tf(gamma,[tau 1]) * tf(wn^2,[1 2*xi*wn wn^2]);

% Add unmodeled dynamics and set SampleStateDim to 5 to get representative
% sample values of the uncertain model P
delta = ultidyn('delta',[1 1],'SampleStateDim',5,'Bound',1);
W = makeweight(0.1,20,10);
P = P * (1+W*delta)

P =

  Uncertain continuous-time state-space model with 1 outputs, 1 inputs, 4 states.
  The model uncertainty consists of the following blocks:
    delta: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences
    gamma: Uncertain real, nominal = 2, variability = [-30,30]%, 1 occurrences
    tau: Uncertain real, nominal = 1, variability = [-30,30]%, 1 occurrences

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

Набор переходных характеристик иллюстрирует изменчивость объекта.

step(P,4)
title('Sampled Step Responses of Uncertain System')

Неопределенная модель объекта управления P содержит 3 неопределенных элемента. Для целей системы управления часто желательно упростить эту модель неопределенности, сохраняя при этом ее общую изменчивость. Это другое использование команды ucover.

Как использовать ucover в этом контексте сначала отобразите неопределенную модель P в массив моделей LTI с использованием usample. Эта команда отсчитывает неопределенные элементы в неопределенной системе и возвращает соответствующие модели LTI, каждая модель представляет одно возможное поведение неопределенной системы. В этом примере выборка P в 60 точках (генератор случайных чисел посеян для повторяемости):

rng(0,'twister');
Parray = usample(P,60);

Далее используйте ucover чтобы охватить все виды поведения в Parray простой моделью неопределенности Usys. Выберите номинальное значение P в качестве центра крышки и использовать фильтр 2-го порядка для моделирования частотного распределения немодулированной динамики.

orderWt = 2;
Parrayg = frd(Parray,logspace(-3,3,60));
[Usys,Info] = ucover(Parrayg,P.NominalValue,orderWt,'InputMult');

График величины Bode показов, как фильтр величины (красным цветом) «покрывает» относительную изменчивость частотной характеристики объекта (синим цветом).

Wt = Info.W1;
bodemag((P.NominalValue-Parray)/P.NominalValue,'b--',Wt,'r')
title('Relative Gaps (blue) vs. Shaping Filter Magnitude (red)')

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

Корректировка взвешивания неопределенности

В этом третьем примере мы начинаем с 40 частотных характеристик 2-входной, 2-выходной системы. Эти данные были собраны с помощью анализатора частоты при различных рабочих условиях. Номинальная модель с двумя состояниями подобрана к наиболее типичному отклику:

A = [-5 10;-10 -5];
B = [1 0;0 1];
C = [1 10;-10 1];
D = 0;
Pnom = ss(A,B,C,D);

Данные частотной характеристики загружаются в массив 40 на 1 моделей FRD:

load ucover_demo
size(Pdata)

40x1 array of FRD models.
Each model has 2 outputs, 2 inputs, and 120 frequency points.

Постройте график этих данных и наложите номинальную модель.

bode(Pdata,'b--',Pnom,'r',{.1,1e3}), grid
legend('Frequency response data','Nominal model','Location','NorthEast')

Поскольку изменчивость отклика скромна, попробуйте смоделировать это семейство частотных характеристик с помощью аддитивной модели неопределенности вида

P = Pnom + w * Delta

где Delta представляет собой ultidyn 2 на 2 объект, представляющий немоделированную динамику и w - скаляр весовая функция, отражающая частотное распределение неопределенности (изменчивость в Pdata).

Начните с фильтра первого порядка w и сравните величину w с минимальной степенью неопределенности, необходимой на каждой частоте:

[P1,InfoS1] = ucover(Pdata,Pnom,1,'Additive');
w = InfoS1.W1;
bodemag(w,'r',InfoS1.W1opt,'g',{1e-1 1e3})
title('Scalar Additive Uncertainty Model')
legend('First-order w','Min. uncertainty amount','Location','SouthWest')

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

[P3,InfoS3] = ucover(Pnom,InfoS1,3,'Additive');
w = InfoS3.W1;
bodemag(w,'r',InfoS3.W1opt,'g',{1e-1 1e3})
title('Scalar Additive Uncertainty Model')
legend('Third-order w','Min. uncertainty amount','Location','SouthWest')

Величина w теперь близко совпадает с минимальной суммой неопределенности. Среди аддитивных моделей неопределенности P3 обеспечивает плотное покрытие поведения в Pdata. Обратите внимание, что P3 имеет в общей сложности 8 состояний (2 от номинальной части и 6 от w).

P3

P3 =

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

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

Можно уточнить эту аддитивную модель неопределенности, используя некалярные функции взвешивания неопределенности, для примера

P = Pnom + W1*Delta*W2

где W1 и W2 представляют собой диагональные фильтры 2 на 2. В этом примере ограничьте использование W2=1 и позволяет обоим диагональным элементам W1 быть третьим порядком.

[PM,InfoM] = ucover(Pdata,Pnom,[3;3],[],'Additive');

Сравните две записи W1 с минимальной суммой неопределенности, рассчитанной ранее. Обратите внимание, что на всех частотах один из диагональных элементов W1 имеет величину, намного меньшую, чем скалярный фильтр w. Это предполагает, что диагонально взвешенная модель неопределенности дает менее консервативное покрытие семейства частотной характеристики.

bodemag(InfoS1.W1opt,'g*',...
    InfoM.W1opt(1,1),'r--',InfoM.W1(1,1),'r',...
    InfoM.W1opt(2,2),'b--',InfoM.W1(2,2),'b',{1e-1 1e3});
title('Diagonal Additive Uncertainty Model')
legend('Scalar Optimal Weight',...
    'W1(1,1), pointwise optimal',...
    'W1(1,1), 3rd-order fit',...
    'W1(2,2), pointwise optimal',...
    'W1(2,2), 3rd-order fit',...
    'Location','SouthWest')

Степень консервативности одного покрытия над другим может быть частично определена количественно с учетом двух частотно-зависимых величин:

Fd2s = norm(inv(W1)*w) ,   Fs2d = norm(W1/w)

Эти величины измеряются тем, насколько одна модель неопределенности должна быть масштабирована, чтобы покрыть другую. Для примера, модель неопределенности Pnom + W1*Delta необходимо увеличить на множитель Fd2s включить все модели, представленные неопределенной моделью Pnom + w*Delta.

График Fd2s и Fs2d как функцию частоты.

Fd2s = fnorm(InfoS1.W1opt*inv(InfoM.W1opt));
Fs2d = fnorm(InfoM.W1opt*inv(InfoS1.W1opt));
semilogx(fnorm(Fd2s),'b',fnorm(Fs2d),'r'), grid
axis([0.1 1000 0.5 2.6])
xlabel('Frequency (rad/s)'), ylabel('Magnitude')
title('Scale factors relating different covers')
legend('Diagonal to Scalar factor',...
    'Scalar to Diagonal factor','Location','SouthWest');

Этот график показывает, что:

Это указывает, что Pnom+W1*Delta модель обеспечивает менее консервативное покрытие данных частотной характеристики в Pdata.