Этот пример показывает, как спроектировать и настроить контроллер с запланированным усилением для химического реактора, переходящего от низкой скорости преобразования к высокой. Для справки см. Seborg, D.E. et al., «Process Dynamics and Control», 2nd Ed., 2004, Wiley, pp. 34-36.
Рассматриваемый здесь способ представляет собой непрерывный реактор бака перемешивания (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, не используя 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, чтобы симулировать нелинейную характеристику контроллера, запланированного по усилению.
setBlockParam
| slTuner
| tunableSurface