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

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

Процесс, рассмотренный здесь, является непрерывным реактором смесителя (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);

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

Чтобы достигнуть этих требований, мы используем контроллер PI во внешнем цикле и ведущий компенсатор во внутреннем цикле. Из-за медленного уровня выборки, ведущий компенсатор необходим, чтобы соответственно стабилизировать химическую реакцию при средней концентрации Cr = 5.28 kmol/m^3/min. Поскольку движущие силы реакции отличаются существенно с концентрацией, мы далее планируем усиления контроллера как функцию концентрации. Это моделируется в 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.21, Hard = 0.99882, Iterations = 205

Получившийся проект удовлетворяет трудное ограничение (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.21, Hard = 0.99988, Iterations = 211

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