Описание объекта

Объект в нашем примере состоит из двух баков с водой в каскаде, как схематически показано на Фигуру 1. Верхний бак (бак 1) питается горячей и холодной водой через управляемые компьютером клапаны. Нижний бак (бак 2) подается водой из выхода в нижней части бака 1. Переполнение поддерживает постоянный уровень в баке 2. Поток смещения холодной воды также подает бак 2 и позволяет бакам иметь различные установившиеся температуры.

Наша цель проекта состоит в том, чтобы контролировать температуры обоих баков 1 и 2. Контроллер имеет доступ к ссылочным командам и измерениям температуры.

Фигура 1: Принципиальная схема системы с двумя резервуарами

Переменные баки

Дадим объект переменным следующие обозначения:

Для удобства зададим систему нормированных модулей следующим образом:

  Variable      Unit Name      0 means:             1 means:
  --------      ---------      --------             --------
  temperature   tunit          cold water temp.     hot water temp.
  height        hunit          tank empty           tank full
  flow          funit          zero input flow      max. input flow

Используя вышеуказанные модули, это параметры объекта управления:

A1 = 0.0256;	% Area of tank 1 (hunits^2)
A2 = 0.0477;	% Area of tank 2 (hunits^2)
h2 = 0.241;	    % Height of tank 2, fixed by overflow (hunits)
fb = 3.28e-5;   % Bias stream flow (hunits^3/sec)
fs = 0.00028;	% Flow scaling (hunits^3/sec/funit)
th = 1.0;	    % Hot water supply temp (tunits)
tc = 0.0;	    % Cold water supply temp (tunits)
tb = tc;	    % Cold bias stream temp (tunits)
alpha = 4876;   % Constant for flow/height relation (hunits/funits)
beta = 0.59;    % Constant for flow/height relation (hunits)

Переменная fs является масштабирующим коэффициент потока, который преобразует вход (0 в 1 фуниты) в поток в hunits ^ 3/second. Константы альфа и бета описывают отношение расход/высота для бака 1:

h1 = альфа * f1-бета.

Номинальные модели бака

Мы можем получить номинальные модели бака путем линеаризации вокруг следующей рабочей точки (все нормированные значения):

h1ss = 0.75;                            % Water level for tank 1
t1ss = 0.75;                            % Temperature of tank 1
f1ss = (h1ss+beta)/alpha;               % Flow tank 1 -> tank 2
fss = [th,tc;1,1]\[t1ss*f1ss;f1ss];
fhss = fss(1);                          % Hot flow
fcss = fss(2);                          % Cold flow
t2ss = (f1ss*t1ss + fb*tb)/(f1ss + fb); % Temperature of tank 2

Номинальная модель для бака 1 имеет входы [fh; fc] и выходы [h1; t1]:

A = [ -1/(A1*alpha),          0;
      (beta*t1ss)/(A1*h1ss),  -(h1ss+beta)/(alpha*A1*h1ss)];

B = fs*[ 1/(A1*alpha),   1/(A1*alpha);
         th/A1,          tc/A1];

C = [ alpha,             0;
      -alpha*t1ss/h1ss,  1/h1ss];

D = zeros(2,2);
tank1nom = ss(A,B,C,D,'InputName',{'fh','fc'},'OutputName',{'h1','t1'});

clf
step(tank1nom), title('Step responses of Tank 1')

Фигура 2: Переходные характеристики бака 1.

Номинальная модель для бака 2 имеет входы [|h1|;|t1|] и выход t2:

A = -(h1ss + beta + alpha*fb)/(A2*h2*alpha);
B = [ (t2ss+t1ss)/(alpha*A2*h2),  (h1ss + beta)/(alpha*A2*h2) ];
C = 1;
D = zeros(1,2);

tank2nom = ss(A,B,C,D,'InputName',{'h1','t1'},'OutputName','t2');

step(tank2nom), title('Step responses of Tank 2')

Фигура 3: Переходные характеристики бака 2.

Модели привода

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

Мы первоначально настраивали модель привода с одним входом (команда сигнал) и двумя выходами (приведенный в действие сигнал и его производная). Мы будем использовать производный выход в ограничении скорости приведения в действие при разработке закона управления.

act_BW = 20;		% Actuator bandwidth (rad/sec)
actuator = [ tf(act_BW,[1 act_BW]); tf([act_BW 0],[1 act_BW]) ];
actuator.OutputName = {'Flow','Flow rate'};

bodemag(actuator)
title('Valve actuator dynamics')

hot_act = actuator;
set(hot_act,'InputName','fhc','OutputName',{'fh','fh_rate'});
cold_act =actuator;
set(cold_act,'InputName','fcc','OutputName',{'fc','fc_rate'});

Фигура 4: Динамика привода клапана.

Сглаживающие фильтры

Все измеренные сигналы фильтруются фильтрами Баттерворта четвертого порядка, каждый с частотой отключения 2,25 Гц.

fbw = 2.25;		% Anti-aliasing filter cut-off (Hz)
filter = mkfilter(fbw,4,'Butterw');
h1F = filter;
t1F = filter;
t2F = filter;

Неопределенность в динамике модели

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

Для примера разомкнутого контура эксперименты указывают на значительное количество динамической неопределенности в t1 ответ. Это связано, в основном, с перемешиванием и потерями тепла. Мы можем смоделировать его как мультипликативную (относительную) ошибку модели, Delta2 в t1 выход. Точно так же мы можем добавить мультипликативные ошибки модели Delta1 и Delta3 к h1 и t2 выходы, показанные на фигуре 5.

Фигура 5: Схематическое представление возмущенной линейной двухбаковой системы.

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

Эти данные предполагают следующие варианты для частотно-зависимых границ ошибки моделирования.

Wh1 = 0.01+tf([0.5,0],[0.25,1]);
Wt1 = 0.1+tf([20*h1ss,0],[0.2,1]);
Wt2 = 0.1+tf([100,0],[1,21]);

clf
bodemag(Wh1,Wt1,Wt2), title('Relative bounds on modeling errors')
legend('h1 dynamics','t1 dynamics','t2 dynamics','Location','NorthWest')

Фигура 6. Относительные ограничения ошибок моделирования.

Теперь мы готовы создать неопределенные модели бака, которые захватывают ошибки моделирования, обсужденные выше.

% Normalized error dynamics
delta1 = ultidyn('delta1',[1 1]);
delta2 = ultidyn('delta2',[1 1]);
delta3 = ultidyn('delta3',[1 1]);

% Frequency-dependent variability in h1, t1, t2 dynamics
varh1 = 1+delta1*Wh1;
vart1 = 1+delta2*Wt1;
vart2 = 1+delta3*Wt2;

% Add variability to nominal models
tank1u = append(varh1,vart1)*tank1nom;
tank2u = vart2*tank2nom;

tank1and2u = [0 1; tank2u]*tank1u;

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

step(tank1u,1000), title('Variability in responses due to modeling errors (Tank 1)')

Фигура 7: Изменчивость откликов из-за ошибок моделирования (Бак 1).

Настройка проектирования контроллера

Теперь давайте рассмотрим задачу системы управления. Мы заинтересованы в отслеживании команд уставки для t1 и t2. Чтобы использовать преимущества алгоритмов разработки H-бесконечности, мы должны сформулировать проект как задачу минимизации усиления с обратной связью. Для этого мы выбираем функции взвешивания, которые захватывают характеристики нарушения порядка и требования к эффективности, чтобы помочь нормализовать соответствующие частотно-зависимые ограничения усиления.

Вот подходящая взвешенная передаточная функция без разомкнутого контура для задачи с двумя резервуарами:

Фигура 8: Система управления соединение для двухбаковой системы.

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

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

Wh1noise = zpk(0.01);  % h1 noise weight
Wt1noise = zpk(0.03);  % t1 noise weight
Wt2noise = zpk(0.03);  % t2 noise weight

Веса ошибок штрафуют ошибки отслеживания уставки на t1 и t2. Мы выберем lowpass первого порядка для этих весов. Мы используем более высокий вес (лучшее отслеживание) для t1 потому что физические факторы заставляют нас верить в это t1 легче контролировать, чем t2.

Wt1perf = tf(100,[400,1]);	% t1 tracking error weight
Wt2perf = tf(50,[800,1]);	% t2 tracking error weight

clf
bodemag(Wt1perf,Wt2perf)
title('Frequency-dependent penalty on setpoint tracking errors')
legend('t1','t2')

Фигура 9: Частотно-зависимый штраф от ошибок отслеживания уставки.

Ссылки (уставки) отражают частотное содержимое таких команд. Поскольку большая часть воды, поступающей в бак 2, поступает из бака 1, изменяется в t2 преобладают изменения в t1. Также t2 обычно командуется значение, близкое к t1. Поэтому имеет больше смысла использовать взвешивание уставки, выраженное в терминах t1 и t2-t1:

  t1cmd = Wt1cmd * w1

  t2cmd = Wt1cmd * w1 + Wtdiffcmd * w2

где w1, w2 - входы белого шума. Адекватный выбор веса:

Wt1cmd = zpk(0.1);               % t1 input command weight
Wtdiffcmd = zpk(0.01);           % t2 - t1  input command weight

Наконец, мы хотели бы наказать и амплитуду, и скорость привода. Мы делаем это путем взвешивания fhc (и fcc) с функцией, которая накатывается на высоких частотах. Кроме того, мы можем создать модель привода с fh и d'fh|/dt как выходы и взвешивать каждый выход отдельно с постоянными весами. Этот подход имеет преимущество уменьшения количества состояний в взвешенной модели разомкнутого контура.

Whact =  zpk(0.01);  % Hot actuator penalty
Wcact =  zpk(0.01);  % Cold actuator penalty

Whrate = zpk(50);    % Hot actuator rate penalty
Wcrate = zpk(50);    % Cold actuator rate penalty

Создание взвешенной модели разомкнутого контура

Теперь, когда мы смоделировали все компоненты объекта и выбрали наши веса проекта, мы будем использовать connect функция для построения неопределенной модели взвешенной модели без разомкнутого контура, показанной на фигуре 8.

inputs = {'t1cmd', 'tdiffcmd', 't1noise', 't2noise', 'fhc', 'fcc'};
outputs = {'y_Wt1perf', 'y_Wt2perf', 'y_Whact', 'y_Wcact', ...
             'y_Whrate', 'y_Wcrate', 'y_Wt1cmd', 'y_t1diffcmd', ...
                                           'y_t1Fn', 'y_t2Fn'};

hot_act.InputName = 'fhc'; hot_act.OutputName = {'fh' 'fh_rate'};
cold_act.InputName = 'fcc'; cold_act.OutputName = {'fc' 'fc_rate'};

tank1and2u.InputName = {'fh','fc'};
tank1and2u.OutputName = {'t1','t2'};

t1F.InputName = 't1'; t1F.OutputName = 'y_t1F';
t2F.InputName = 't2'; t2F.OutputName = 'y_t2F';

Wt1cmd.InputName = 't1cmd'; Wt1cmd.OutputName = 'y_Wt1cmd';
Wtdiffcmd.InputName = 'tdiffcmd'; Wtdiffcmd.OutputName = 'y_Wtdiffcmd';

Whact.InputName = 'fh'; Whact.OutputName = 'y_Whact';
Wcact.InputName = 'fc'; Wcact.OutputName = 'y_Wcact';

Whrate.InputName = 'fh_rate'; Whrate.OutputName = 'y_Whrate';
Wcrate.InputName = 'fc_rate'; Wcrate.OutputName = 'y_Wcrate';

Wt1perf.InputName = 'u_Wt1perf'; Wt1perf.OutputName = 'y_Wt1perf';
Wt2perf.InputName = 'u_Wt2perf'; Wt2perf.OutputName = 'y_Wt2perf';

Wt1noise.InputName = 't1noise'; Wt1noise.OutputName = 'y_Wt1noise';
Wt2noise.InputName = 't2noise'; Wt2noise.OutputName = 'y_Wt2noise';

sum1 = sumblk('y_t1diffcmd = y_Wt1cmd + y_Wtdiffcmd');
sum2 = sumblk('y_t1Fn = y_t1F + y_Wt1noise');
sum3 = sumblk('y_t2Fn = y_t2F + y_Wt2noise');
sum4 = sumblk('u_Wt1perf = y_Wt1cmd - t1');
sum5 = sumblk('u_Wt2perf = y_Wtdiffcmd + y_Wt1cmd - t2');

% This produces the uncertain state-space model
P = connect(tank1and2u,hot_act,cold_act,t1F,t2F,Wt1cmd,Wtdiffcmd,Whact, ...
                Wcact,Whrate,Wcrate,Wt1perf,Wt2perf,Wt1noise,Wt2noise, ...
                   sum1,sum2,sum3,sum4,sum5,inputs,outputs);

disp('Weighted open-loop model: ')
P

Weighted open-loop model: 

P =

  Uncertain continuous-time state-space model with 10 outputs, 6 inputs, 18 states.
  The model uncertainty consists of the following blocks:
    delta1: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences
    delta2: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences
    delta3: 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.

Проектирование контроллера H-бесконечности

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

nmeas = 4;		% Number of measurements
nctrls = 2;		% Number of controls
[k0,g0,gamma0] = hinfsyn(P.NominalValue,nmeas,nctrls);
gamma0

gamma0 =

    0.9016

Наименьший достижимый коэффициент усиления в системе с обратной связью составляет около 0,9, что показывает нам, что наши спецификации эффективности отслеживания частотного диапазона выполняются контроллером k0. Симуляция этого проекта во временном интервале является разумным способом проверить, что мы правильно установили веса эффективности. Во-первых, мы создадим модель с обратной связью, отображающую входные сигналы [t1ref; t2ref; t1noise; t2noise] к выходным сигналам [h1; t1; t2; fhc; fcc]:

inputs = {'t1ref', 't2ref', 't1noise', 't2noise', 'fhc', 'fcc'};
outputs = {'y_tank1', 'y_tank2', 'fhc', 'fcc', 'y_t1ref', 'y_t2ref', ...
                'y_t1Fn', 'y_t2Fn'};

hot_act(1).InputName = 'fhc'; hot_act(1).OutputName = 'y_hot_act';
cold_act(1).InputName = 'fcc'; cold_act(1).OutputName = 'y_cold_act';

tank1nom.InputName = [hot_act(1).OutputName cold_act(1).OutputName];
tank1nom.OutputName = 'y_tank1';
tank2nom.InputName = tank1nom.OutputName;
tank2nom.OutputName = 'y_tank2';

t1F.InputName = tank1nom.OutputName(2); t1F.OutputName = 'y_t1F';
t2F.InputName = tank2nom.OutputName; t2F.OutputName = 'y_t2F';

I_tref = zpk(eye(2));
I_tref.InputName = {'t1ref', 't2ref'}; I_tref.OutputName = {'y_t1ref', 'y_t2ref'};

sum1 = sumblk('y_t1Fn = y_t1F + t1noise');
sum2 = sumblk('y_t2Fn = y_t2F + t2noise');

simlft = connect(tank1nom,tank2nom,hot_act(1),cold_act(1),t1F,t2F,I_tref,sum1,sum2,inputs,outputs);

% Close the loop with the H-infinity controller |k0|
sim_k0 = lft(simlft,k0);
sim_k0.InputName = {'t1ref'; 't2ref'; 't1noise'; 't2noise'};
sim_k0.OutputName = {'h1'; 't1'; 't2'; 'fhc'; 'fcc'};

Теперь моделируем реакцию с обратной связью при нарастании уставок на t1 и t2 от 80 секунд до 100 секунд:

time=0:800;
t1ref = (time>=80 & time<100).*(time-80)*-0.18/20 + ...
    (time>=100)*-0.18;
t2ref = (time>=80 & time<100).*(time-80)*-0.2/20 + ...
    (time>=100)*-0.2;
t1noise = Wt1noise.k * randn(size(time));
t2noise = Wt2noise.k * randn(size(time));

y = lsim(sim_k0,[t1ref ; t2ref ; t1noise ; t2noise],time);

Далее добавим моделируемые выходы к их устойчивым состояниям значениям и постройм график откликов:

h1 = h1ss+y(:,1);
t1 = t1ss+y(:,2);
t2 = t2ss+y(:,3);
fhc = fhss/fs+y(:,4); % Note scaling to actuator
fcc = fcss/fs+y(:,5); % Limits (0<= fhc <= 1) etc.

В этом коде мы строим графики выходов t1 и t2, а также высота h1 бака 1:

plot(time,h1,'--',time,t1,'-',time,t2,'-.');
xlabel('Time (sec)')
ylabel('Measurements')
title('Step Response of H-infinity Controller k0')
legend('h1','t1','t2');
grid

Фигура 10: Переходная характеристика контроллера k0 H-бесконечности.

Далее строим графики команд на горячий и холодный приводы.

plot(time,fhc,'-',time,fcc,'-.');
xlabel('Time: seconds')
ylabel('Actuators')
title('Actuator Commands for H-infinity Controller k0')
legend('fhc','fcc');
grid

Фигура 11: Команды привода для контроллера H-бесконечности k0.

Робастность контроллера H-бесконечности

Контроллер H-бесконечности k0 предназначен для номинальных моделей бака. Давайте рассмотрим, насколько хорошо ее тарифы для возмущенной модели в пределах неопределенности модели. Мы можем сравнить номинальную эффективность с обратной связью gamma0 с наихудшей эффективностью по сравнению с набором неопределенностей модели. (для получения дополнительной информации см. неопределенность в динамике модели»).

clpk0 = lft(P,k0);

% Compute and plot worst-case gain
wcsigmaplot(clpk0,{1e-4,1e2})
ylim([-20 10])

Фигура 12: Анализ эффективности для контроллера k0.

Худшая эффективность замкнутого контура значительно хуже номинальной эффективности, которая говорит нам, что контроллер H-бесконечности k0 не является устойчивым к ошибкам моделирования.

Синтез контроллера Mu

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

[kmu,bnd] = musyn(P,nmeas,nctrls);

D-K ITERATION SUMMARY:
-----------------------------------------------------------------
                       Robust performance               Fit order
-----------------------------------------------------------------
  Iter         K Step       Peak MU       D Fit             D
    1           8.814        2.862        2.888             4
    2           2.497        2.092        2.112            10
    3           1.351        1.329        1.337            10
    4           1.201        1.198        1.209            16
    5           1.191         1.19        1.196            16
    6           1.192         1.19        1.196            20

Best achieved robust performance: 1.19

Как и прежде, мы можем симулировать отклики с обратной связью с контроллером kmu

sim_kmu = lft(simlft,kmu);
y = lsim(sim_kmu,[t1ref;t2ref;t1noise;t2noise],time);
h1 = h1ss+y(:,1);
t1 = t1ss+y(:,2);
t2 = t2ss+y(:,3);
fhc = fhss/fs+y(:,4); % Note scaling to actuator
fcc = fcss/fs+y(:,5); % Limits (0<= fhc <= 1) etc.

% Plot |t1| and |t2| as well as the height |h1| of tank 1
plot(time,h1,'--',time,t1,'-',time,t2,'-.');
xlabel('Time: seconds')
ylabel('Measurements')
title('Step Response of mu Controller kmu')
legend('h1','t1','t2');
grid

Фигура 13: Переходная характеристика mu контроллера kmu.

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

% Worst-case performance for kmu
clpmu = lft(P,kmu);
wcsigmaplot(clpmu,{1e-4,1e2})
ylim([-20 10])

Фигура 14: Анализ эффективности для контроллера kmu.

Можно использовать wcgain непосредственно вычислить коэффициент усиления в худшем случае на частоте (наихудший случай пикового усиления или наихудшая норма H-бесконечности). Можно также вычислить его чувствительность к каждому неопределенному элементу. Результаты показывают, что пик усиления в худшем случае наиболее чувствителен к изменениям в области значений delta2.

opt = wcOptions('Sensitivity','on');
[wcg,wcu,wcinfo] = wcgain(clpmu,opt);
wcg

wcg = 

  struct with fields:

           LowerBound: 1.3143
           UpperBound: 1.3171
    CriticalFrequency: 0

wcinfo.Sensitivity

ans = 

  struct with fields:

    delta1: 0
    delta2: 59
    delta3: 10