Реактор с непрерывным перемешиванием Бака

Рассматриваемый здесь способ представляет собой непрерывный реактор бака перемешивания (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';

Блок Предел на контроллер выходе задает, что температура хладагента 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);

Планируемый по усилению контроллер

Для достижения этих требований мы используем ПИ-контроллер во внешнем контуре и свинцовый компенсатор во внутреннем цикле. Из-за медленной частоты дискретизации, свинцовый компенсатор необходим для адекватной стабилизации химической реакции при средней концентрации Cr = 5,28 кмоль/м ^ 3/мин. Поскольку динамика реакции существенно изменяется в зависимости от концентрации, мы далее планируем коэффициент усиления контроллера как функцию концентрации. Это смоделировано в Simulink с использованием блоков Интерполяционной таблицы, как показано фигуры и 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);

Настройка контроллера

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