Документация

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

Объект в нашем примере состоит из двух баков для воды в каскаде как показано схематично в рисунке 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; 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: Схематическое представление встревоженной, линейной системы 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. Мы выберем фильтры нижних частот первого порядка для этих весов. Мы используем более высокий вес (лучше отслеживающий) для 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
wcsigma(clpk0)
axis([1e-4 100 -20 10])

Рисунок 12: анализ Производительности для контроллера k0.

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

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

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

% Options for D-K iterations
frad = 2*pi*logspace(-5,1,30);
opt = dksynOptions('FrequencyVector',frad,'NumberofAutoIterations',4);

% Perform mu synthesis
[kmu,clpmu,bnd] = dksyn(P,nmeas,nctrls,opt);
bnd

bnd =

    1.2055

Как прежде, мы можем моделировать ответы с обратной связью с контроллером 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
wcsigma(clpmu)
axis([1e-4 100 -20 10])

Рисунок 14: анализ Производительности для контроллера kmu.

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

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

wcg = 

  struct with fields:

           LowerBound: 1.3250
           UpperBound: 1.3278
    CriticalFrequency: 8.3879e-04

wcinfo.Sensitivity

ans = 

  struct with fields:

    delta1: 0
    delta2: 52
    delta3: 5