Описание установки

Установка в нашем примере состоит из двух резервуаров для воды в каскаде, как схематически показано на фиг.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. Мы выберем фильтры нижних частот первого порядка для этих весов. Мы используем более высокий вес (лучшее отслеживание) для 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: Пошаговый отклик контроллера H-бесконечности k0.

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

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