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

Процесс, рассмотренный здесь, является непрерывным реактором смесителя (CSTR) во время перехода от низко до высокой скорости преобразования (высоко к низкой остаточной концентрации). Поскольку химическая реакция является экзотермической (производит тепло), реакторная температура должна контролироваться, чтобы предотвратить тепловой пробой. Задача управления осложнена тем, что движущие силы процесса нелинейны и переход от устойчивого до нестабильного и назад в стабильное когда скорость преобразования увеличивается. Реакторные движущие силы моделируются в Simulink. Управляемые переменные являются остаточной концентрацией Cr и реакторный температурный Tr, и переменной, которой управляют, является температурный Tc из хладагента, циркулирующего в конверте охлаждения реактора.

open_system('rct_CSTR_OL')

Мы хотим перейти от остаточной концентрации 8.57 kmol/m^3 первоначально вниз к 2 kmol/m^3. Изучать, как движущие силы процесса развивают с остаточной концентрацией Cr, найдите условия равновесия для пяти значений Cr между 8,57 и 2 и линеаризуют динамику процесса вокруг каждого равновесия. Регистрируйте реактор и температуры хладагента в каждой точке равновесия.

CrEQ = linspace(8.57,2,5)';  % concentrations
TrEQ = zeros(5,1);           % reactor temperatures
TcEQ = zeros(5,1);           % coolant temperatures

% Specify trim conditions
opspec = operspec('rct_CSTR_OL',5);
for k=1:5
   % Set desired residual concentration
   opspec(k).Outputs(1).y = CrEQ(k);
   opspec(k).Outputs(1).Known = true;
end

% Compute equilibrium condition and log corresponding temperatures
[op,report] = findop('rct_CSTR_OL',opspec,...
   findopOptions('DisplayReport','off'));
for k=1:5
   TrEQ(k) = report(k).Outputs(2).y;
   TcEQ(k) = op(k).Inputs.u;
end

% Linearize process dynamics at trim conditions
G = linearize('rct_CSTR_OL', 'rct_CSTR_OL/CSTR', op);
G.InputName = {'Cf','Tf','Tc'};
G.OutputName = {'Cr','Tr'};

Постройте реактор и температуры хладагента в равновесии в зависимости от концентрации.

subplot(311), plot(CrEQ,'b-*'), grid, title('Residual concentration'), ylabel('CrEQ')
subplot(312), plot(TrEQ,'b-*'), grid, title('Reactor temperature'), ylabel('TrEQ')
subplot(313), plot(TcEQ,'b-*'), grid, title('Coolant temperature'), ylabel('TcEQ')

Стратегия регулирования без обратной связи состоит из следующих профиль температуры хладагента выше к гладко переходу между Cr=8.57 и Cr=2 равновесия. Однако эта стратегия обречена тем, что реакция нестабильна в среднем и должна правильно охлаждаться, чтобы избежать теплового пробоя. Это подтверждено путем осмотра полюсов линеаризовавших моделей для этих пяти точек равновесия, рассмотренных выше (три из этих пяти моделей, нестабильны).

pole(G)

ans(:,:,1) =

  -0.5225 + 0.0000i
  -0.8952 + 0.0000i


ans(:,:,2) =

   0.1733 + 0.0000i
  -0.8866 + 0.0000i


ans(:,:,3) =

   0.5114 + 0.0000i
  -0.8229 + 0.0000i


ans(:,:,4) =

   0.0453 + 0.0000i
  -0.4991 + 0.0000i


ans(:,:,5) =

  -1.1077 + 1.0901i
  -1.1077 - 1.0901i

Диаграмма Боде дальнейшие подсветки значительные изменения динамики объекта при переходе от Cr=8.57 к Cr=2.

clf, bode(G(:,'Tc'),{0.01,10})

Стратегия управления с обратной связью

Чтобы предотвратить тепловой пробой при сползании вниз остаточной концентрации, используйте управление с обратной связью, чтобы скорректировать температуру хладагента Tc на основе измерений остаточной концентрации Cr и реакторный температурный Tr. Для этого приложения мы используем архитектуру каскадного регулирования, где внутренний цикл регулирует реакторную температуру, и внешний контур отслеживает заданное значение концентрации. Обе обратной связи является цифровой с периодом выборки 0,5 минут.

open_system('rct_CSTR')

Целевая концентрация Cref пандусы по сравнению с 8.57 kmol/m^3 в t=10 к 2 kmol/m^3 в t=36 (переход длится 26 минут). Соответствующий профиль Tref поскольку реакторная температура получена путем интерполяции значений равновесия TrEQ от анализа для обрезки. Контроллер вычисляет корректировку температуры хладагента dTc относительно начального значения равновесия TcEQ(1)=297.98 для Cr=8.57. Обратите внимание на то, что модель настраивается так, чтобы первоначально, выход TrSP из "Concentration controller" блок совпадает с реакторной температурой, корректировка dTc нуль и температура хладагента Tc в его значении равновесия TcEQ(1).

clf
t = [0 10:36 45];
C = interp1([0 10 36 45],[8.57 8.57 2 2],t);
subplot(211), plot(t,C), grid, set(gca,'ylim',[0 10])
title('Target residual concentration'), ylabel('Cref')
subplot(212), plot(t,interp1(CrEQ,TrEQ,C))
title('Corresponding reactor temperature at equilibrium'), ylabel('Tref'), grid

Управляйте целями

Используйте TuningGoal объекты получить конструктивные требования. Во-первых, Cr должен следовать за заданными значениями Cref со временем отклика приблизительно в 5 минут.

R1 = TuningGoal.Tracking('Cref','Cr',5);

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

MinDecay = 0.2;
MinDamping = 0.5;
% Constrain closed-loop poles of inner loop with the outer loop open
R2 = TuningGoal.Poles('Tc',MinDecay,MinDamping);
R2.Openings = 'TrSP';

Блок Rate Limit в контроллере выход указывает что температура хладагента Tc не может варьироваться быстрее, чем 10 градусов в минуту. Это - серьезное ограничение на полномочия контроллера, которые, когда проигнорировано, могут привести к низкой производительности или нестабильности. Чтобы принять это ограничение скорости во внимание, наблюдайте тот Cref варьируется на уровне 0.25 kmol/m^3/min. Гарантировать тот Tc не варьируется быстрее, чем 10 степеней/min, усиление от Cref к Tc должен быть меньше 10/0.25=40.

R3 = TuningGoal.Gain('Cref','Tc',40);

Наконец, потребуйте по крайней мере 7 дБ запаса по амплитуде и 45 градусов запаса по фазе во входе Tc объекта.

R4 = TuningGoal.Margins('Tc',7,45);

Запланированный на усиление контроллер

Чтобы достигнуть этих требований, мы используем ПИ-контроллер во внешнем контуре и ведущий компенсатор во внутреннем цикле. Из-за медленной частоты дискретизации, ведущий компенсатор необходим, чтобы соответственно стабилизировать химическую реакцию при средней концентрации Cr = 5.28 kmol/m^3/min. Поскольку движущие силы реакции варьируются существенно с концентрацией, мы далее планируем усиления контроллера в зависимости от концентрации. Это моделируется в Simulink с помощью блоков Интерполяционной таблицы как показано 1 в цифрах и 2.

Рисунок 1: запланированный на усиление ПИ-контроллер для цикла концентрации.

Рисунок 2: запланированный на усиление ведущий компенсатор для температурного цикла.

Настройка этого запланированного на усиление контроллера составляет настройку данных интерполяционной таблицы в области значений значений концентрации. Вместо того, чтобы настраивать отдельные записи интерполяционной таблицы, параметрируйте усиления контроллера Kp,Ki,Kt,a,b как квадратичные полиномы в Cr, например,

Помимо сокращения количества переменных, чтобы настроиться, этот подход гарантирует сглаженные переходы усиления как Cr варьируется. Используя systune, можно автоматически настроить коэффициенты, чтобы удовлетворить требования R1-R4 в этих пяти точках равновесия, вычисленных выше. Это составляет настройку запланированного на усиление контроллера в пяти точках проекта вдоль Cref траектория. Используйте tunableSurface объект параметрировать каждое усиление как квадратичную функцию Cr. "Настраивающаяся сетка" установлена в эти пять концентраций CrEQ и основные функции для квадратичной параметризации. Большинство усилений инициализируется, чтобы быть тождественно нулевым.

TuningGrid = struct('Cr',CrEQ);
ShapeFcn = @(Cr) [Cr , Cr^2];

Kp = tunableSurface('Kp', 0, TuningGrid, ShapeFcn);
Ki = tunableSurface('Ki', -2, TuningGrid, ShapeFcn);
Kt = tunableSurface('Kt', 0, TuningGrid, ShapeFcn);
a = tunableSurface('a', 0, TuningGrid, ShapeFcn);
b = tunableSurface('b', 0, TuningGrid, ShapeFcn);

Контроллер, настраивающийся

Поскольку целевая полоса пропускания в течение десятилетия после частоты Найквиста, легче настроить контроллер непосредственно в дискретной области. Дискретизируйте линеаризовавшую динамику процесса с шагом расчета 0,5 минут. Используйте метод ZOH, чтобы отразить, как цифровой контроллер взаимодействует с объектом непрерывного времени.

Ts = 0.5;
Gd = c2d(G,Ts);

Создайте slTuner интерфейс для настройки квадратичных расписаний усиления, введенных выше. Используйте замену блока, чтобы заменить нелинейную модель объекта управления пятью дискретизированными линейными моделями Gd полученный в проекте указывает CrEQ. Используйте setBlockParam сопоставлять настраиваемое усиление функционирует Kp, Ki, KtAB с блоками Интерполяционной таблицы того же имени.

BlockSubs = struct('Name','rct_CSTR/CSTR','Value',Gd);
ST0 = slTuner('rct_CSTR',{'Kp','Ki','Kt','a','b'},BlockSubs);
ST0.Ts = Ts;  % sample time for tuning

% Register points of interest
ST0.addPoint({'Cref','Cr','Tr','TrSP','Tc'})

% Parameterize look-up table blocks
ST0.setBlockParam('Kp',Kp);
ST0.setBlockParam('Ki',Ki);
ST0.setBlockParam('Kt',Kt);
ST0.setBlockParam('a',a);
ST0.setBlockParam('b',b);

Можно теперь использовать systune настроить коэффициенты контроллера против требований R1-R4. Сделайте требование по запасу устойчивости трудными ограничениями и оптимизируйте остающиеся требования.

ST = systune(ST0,[R1 R2 R3],R4);

Final: Soft = 1.23, Hard = 0.99989, Iterations = 197

Получившийся проект удовлетворяет трудному ограничению (Hard<1) и почти удовлетворяет остающимся требованиям (Soft близко к 1). Чтобы подтвердить этот проект, симулируйте ответы на пандус в концентрации с тем же наклоном как Cref. Каждый график показывает линейные ответы в пяти точках проекта CrEQ.

t = 0:Ts:20;
uC = interp1([0 2 5 20],(-0.25)*[0 0 3 3],t);
subplot(211), lsim(getIOTransfer(ST,'Cref','Cr'),uC)
grid, set(gca,'ylim',[-1.5 0.5]), title('Residual concentration')
subplot(212), lsim(getIOTransfer(ST,'Cref','Tc'),uC)
grid, title('Coolant temperature variation')

Обратите внимание на то, что скорость изменения температуры хладагента остается в физических пределах (10 градусов в минуту или 5 градусов в период расчета).

Валидация контроллера

Смотрите, как каждое усиление меняется в зависимости от Cr во время перехода.

% Access tuned gain schedules
TGS = getBlockParam(ST);

% Plot gain profiles
clf
subplot(321), viewSurf(TGS.Kp), ylabel('Kp')
subplot(322), viewSurf(TGS.Ki), ylabel('Ki')
subplot(323), viewSurf(TGS.Kt), ylabel('Kt')
subplot(324), viewSurf(TGS.a), ylabel('a')
subplot(325), viewSurf(TGS.b), ylabel('b')

Чтобы подтвердить запланированный на усиление контроллер в Simulink, сначала используйте writeBlockValue применять настраивающиеся результаты к модели Simulink. Для каждого блока Lookup Table это оценивает соответствующую квадратичную формулу усиления в табличных точках останова и обновляет табличные данные соответственно.

writeBlockValue(ST)

Затем нажмите на кнопку воспроизведения, чтобы симулировать ответ с настроенными расписаниями усиления. Результаты симуляции появляются в рисунке 3. Запланированный на усиление контроллер успешно управляет реакцией посредством перехода с соответствующим временем отклика и никаким насыщением ограничений скорости (контроллер выход du совпадает с эффективным температурным изменением dTc). Реакторная температура остается близко к ее значению равновесия Tref, указание, что контроллер сохраняет реакцию около равновесия при предотвращении теплового пробоя.

Рисунок 3: Переход с запланированным на усиление каскадным контроллером.

Контроллер, настраивающий MATLAB

В качестве альтернативы можно настроить расписания усиления непосредственно в MATLAB, не используя slTuner интерфейс. Сначала параметрируйте усиления как квадратичные функции Cr как сделано выше.

TuningGrid = struct('Cr',CrEQ);
ShapeFcn = @(Cr) [Cr , Cr^2];

Kp = tunableSurface('Kp', 0, TuningGrid, ShapeFcn);
Ki = tunableSurface('Ki', -2, TuningGrid, ShapeFcn);
Kt = tunableSurface('Kt', 0, TuningGrid, ShapeFcn);
a = tunableSurface('a', 0, TuningGrid, ShapeFcn);
b = tunableSurface('b', 0, TuningGrid, ShapeFcn);

Используйте эти усиления, чтобы создать PI и привести контроллеры.

PI = pid(Kp,Ki,'Ts',Ts,'TimeUnit','min');
PI.u = 'ECr';   PI.y = 'TrSP';

LEAD = Kt * tf([1 -a],[1 -b],Ts,'TimeUnit','min');
LEAD.u = 'ETr';   LEAD.y = 'Tc';

Используйте connect создавать модель с обратной связью полной системы управления в пяти точках проекта. Отметьте контроллеры выход TrSP и Tc как "указывает анализ" так, чтобы циклы могли быть открыты, и запасы устойчивости оценены в этих местоположениях. Модель T0 с обратной связью массив 5 на 1 линейных моделей в зависимости от настраиваемых коэффициентов Kp,Ki,Kt,a,b. Каждая модель дискретна и произвела каждую половину минуты.

Gd.TimeUnit = 'min';
S1 = sumblk('ECr = Cref - Cr');
S2 = sumblk('ETr = TrSP - Tr');
T0 = connect(Gd(:,'Tc'),LEAD,PI,S1,S2,'Cref','Cr',{'TrSP','Tc'});

Наконец, используйте systune чтобы настроить усиление планируют коэффициенты.

T = systune(T0,[R1 R2 R3],R4);

Final: Soft = 1.23, Hard = 0.99987, Iterations = 193

Результат похож на тот, полученный выше. Подтвердите путем графического вывода усилений в зависимости от Cr использование настроенных коэффициентов в T.

clf
subplot(321), viewSurf(setBlockValue(Kp,T)), ylabel('Kp')
subplot(322), viewSurf(setBlockValue(Ki,T)), ylabel('Ki')
subplot(323), viewSurf(setBlockValue(Kt,T)), ylabel('Kt')
subplot(324), viewSurf(setBlockValue(a,T)), ylabel('a')
subplot(325), viewSurf(setBlockValue(b,T)), ylabel('b')

Можно далее подтвердить проект путем симуляции линейных ответов в каждой точке проекта. Однако необходимо возвратиться к Simulink, чтобы симулировать нелинейный ответ запланированного на усиление контроллера.