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

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

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

Рисунок 1: Принципиальная схема системы 2D бака

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

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

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

  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)

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

h1 = alpha*f1-beta.

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

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

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ФК ] и выходные параметры [h1t1 ]:

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: Схематическое представление встревоженной, линейной системы 2D бака.

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

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

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-бесконечность проектируют алгоритмы, мы должны сформулировать проект как проблему минимизации усиления с обратной связью. Для этого мы выбираем функции взвешивания, которые получают характеристики воздействия и требования к производительности помочь нормировать соответствующие зависимые частотой ограничения усиления.

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

Рисунок 8: соединение Системы управления для системы 2D бака.

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

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

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] к выходным сигналам [h1t1 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 не устойчиво к моделированию ошибок.

Контроллер Му синтез

Чтобы исправить это отсутствие робастности, мы будем использовать 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.494        2.088        2.109            10
    3           1.354        1.331        1.339            10
    4           1.201        1.198        1.209            22
    5           1.194        1.191        1.198            20
    6           1.192         1.19        1.194            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.3166
           UpperBound: 1.3193
    CriticalFrequency: 0

wcinfo.Sensitivity

ans = 

  struct with fields:

    delta1: 0
    delta2: 60
    delta3: 10