Реактор с мешалкой непрерывного действия

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

open_system('rct_CSTR_OL')

Мы хотим перейти от остаточной концентрации 8,57 кмоль/м ^ 3 первоначально к 2 кмоль/м ^ 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 кмоль/м ^ 3 при t = 10 на 2 кмоль/м ^ 3 при t = 36 (переход длится 26 мин). Соответствующий профильTref для температуры реактора получают интерполяцией равновесных значений TrEQ из анализа обрезки. Контроллер вычисляет регулировку температуры хладагента dTc относительно начального равновесного значения TcEQ(1)= 297,98 дляCr=8.57. Обратите внимание, что модель настроена так, что первоначально вывод TrSP блока «Регулятор концентрации» соответствует температуре реактора, регулировка 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 кмоль/м ^ 3/мин. Обеспечить, чтобыTc не изменяется быстрее 10 градусов/мин, коэффициент усиления от Cref кому Tc должно быть меньше 10/0,25 = 40.

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

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

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

Контроллер с планированием усиления

Для выполнения этих требований мы используем контроллер PI во внешнем контуре и компенсатор выводов во внутреннем контуре. Из-за низкой скорости отбора проб для адекватной стабилизации химической реакции при средней концентрации необходим свинцовый компенсатор. Cr = 5,28 кмоль/м ^ 3/мин. Поскольку динамика реакции существенно варьирует в зависимости от концентрации, мы дополнительно планируем усиления контроллера как функцию концентрации. Это моделируется в Simulink с использованием блоков таблицы поиска, как показано на рис. 1 и 2.

Рис. 1. PI-контроллер с планированием усиления для контура концентрации.

Рис. 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, Kt, a, b с блоками таблицы подстановки с тем же именем.

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.99987, Iterations = 193

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

writeBlockValue(ST)

Затем нажмите кнопку Play, чтобы смоделировать ответ с настроенными графиками усиления. Результаты моделирования представлены на рис. 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.99991, Iterations = 195

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