Неадиабатический непрерывный реактор смесителя: моделирование файла MATLAB с симуляциями в Simulink®
В этом примере показано, как включать и симулировать модель IDNLGREY в Simulink®. Мы используем систему химической реакции в качестве основания моделирования. Первая часть моделирования и идентификации примера может быть запущена без Simulink.
Моделирование неадиабатического непрерывного реактора смесителя
Довольно общая химическая система, с которой сталкиваются в перерабатывающей промышленности, является Постоянно реактором смесителя (CSTR). Здесь мы изучим покрытое кожухом диабатическое (i.e., неадиабатическая) реактор бака, описанный экстенсивно в книге Бекетт "Динамика Процесса: Моделирование, Анализ и Симуляция", опубликованный Prentice Hall, 1998. Судно принято, чтобы быть отлично смешанным, и происходит одна экзотермическая и необратимая реакция первого порядка,-> B. Принципиальную схему судна и конверта охлаждения окружения показывают в окне графика. Заметьте, что это - эскиз; в действительности поток хладагента, e.g., обычно окружая целый реакторный конверт, и не только нижнюю часть его.
Рисунок 1: Принципиальная схема CSTR.
Модель CSTR требуется для более усовершенствованных подходов управления. Входной поток реагента A питается на постоянном уровне F. После побуждения потоков конечного продукта из судна на том же уровне как реагент A подан в бак (объем V в реакторном баке таким образом сохранен постоянным). Стратегия управления требует, чтобы температурой конверта u_3 (t) управляли для того, чтобы сохранить концентрацию реагента y_1 (t) на желаемом уровне, несмотря на воздействия, являющиеся результатом входной потоковой концентрации канала и температуры (вводит u_1 (t) и u_2 (t), соответственно). Как температура в баке y_2 (t) может варьироваться значительно во время работы реактора, также желательно гарантировать, что эта переменная процесса остается в рамках разумных пределов.
Моделирование CSTR
Система CSTR моделируется с помощью основных принципов учета и энергосбережения. Изменение концентрации реагента в судне на единицу измерения времени d C_A (t)/dt (= d y_1 (t)/dt) может быть смоделировано как:
d C_A(t)
-------- = F/V*(C_Af(t)-C_A(t)) - r(t)
dtгде первый срок описывает изменения концентрации из-за различий между концентрацией реагента во входном потоке и в судне, и второй срок описывает изменения концентрации (скорость реакции), которая происходит из-за химической реакции в судне. Скорость реакции на единичный объем описана законом об уровне Аррениуса:
r(t) = k_0*exp(-E/(R*T(t)))*C_A(t)
который утверждает, что уровень химической реакции увеличивается экспоненциально с абсолютной температурой. k_0 является здесь неизвестной нетепловой константой, E является энергией активации, идеальной газовой константой Р Больцманна и T (t) (= y_2 (t)) температура в реакторе.
Точно так же с помощью энергетического принципа баланса (принимающий постоянный объем в реакторе), изменение температуры на единицу измерения времени d T (t)/dt в реакторе может быть смоделировано как:
d T(t)
------ = F/V(T_f(t)-T(t)) - (H/c_p*rho)*r(t) - (U*A)/(c_p*rho*V)*(T(t)-T_j(t))
dtгде первые и третьи сроки описывают изменения из-за этого, потоковая температура канала T_f (t) и температура хладагента конверта T_j (t) отличаются от реакторной температуры, соответственно. Второй срок является влиянием на реакторную температуру, вызванную химической реакцией в судне. В этом уравнении H является теплом параметра реакции, c_p термин теплоемкости, ро термин плотности, U полный коэффициент теплопередачи и область для теплообмена (область хладагента/судна).
Соедините, CSTR имеет три входных сигнала:
u_1(t) = C_Af(t) Concentration of A in inlet feed stream [kgmol/m^3]. u_2(t) = T_f(t) Inlet feed stream temperature [K]. u_3(t) = T_j(t) Jacket coolant temperature [K].
и два выходных сигнала:
y_1(t) = C_A(t) Concentration of A in reactor tank [kgmol/m^3]. y_2(t) = T(t) Reactor temperature [K].
которые также являются естественными состояниями модели, i.e., y_1 (t) = x_1 (t) и y_2 (t) = x_2 (t).
После смешания некоторых исходных параметров мы заканчиваем с восемью различными параметрами модели:
F Volumetric flow rate (volume/time) [m^3/h]. Fixed.
V Volume in reactor [m^3]. Fixed.
k_0 Pre-exponential nonthermal factor [1/h]. Free.
E Activation energy [kcal/kgmol]. Free.
R Boltzmann's gas constant [kcal/(kgmol*K)]. Fixed.
H Heat of reaction [kcal/kgmol]. Fixed.
HD = c_p*rho Heat capacity times density [kcal/(m^3*K)]. Free.
HA = U*A Overall heat transfer coefficient times tank
area [kcal/(K*h)] Free.Четыре из параметров здесь заданы, чтобы быть свободными (i.e., чтобы быть оцененным). На практике можно было, вероятно, также определить предэкспоненциальный нетепловой фактор k_0 и энергия активации E из лабораторных экспериментов. Это затем упростило бы идентификационную задачу только рассмотреть два неизвестные: теплоемкость c_p и полный коэффициент теплопередачи U (которые выведены из HD и HA, соответственно, как, ро и A известны).
С вышеупомянутым введенным обозначением следующее представление пространства состояний получено для CSTR.
d x_1(t)
-------- = F/V*(u_1(t)-x_1(t)) - k_0*exp(-E/(R*x_2(t)))*x_1(t)
dt d x_2(t)
-------- = F/V(u_2(t)-x_2(t)) - (H/HD)*k_0*exp(-E/(R*x_2(t)))*x_1(t)
dt - (HA/(HD*V))*(x_2(t)-u_3(t)) y_1(t) = x_1(t)
y_2(t) = x_2(t)Создание неадиабатического непрерывного реактора смесителя объект IDNLGREY
Эта информация вводится в файл MATLAB, названный cstr_m.m со следующим содержимым.
function [dx, y] = cstr_m(t, x, u, F, V, k_0, E, R, H, HD, HA, varargin) %CSTR_M A non-adiabatic Continuous Stirred Tank Reactor (CSTR).
% Output equations.
y = [x(1); ... % Concentration of substance A in the reactor.
x(2) ... % Reactor temperature.
]; % State equations.
dx = [F/V*(u(1)-x(1))-k_0*exp(-E/(R*x(2)))*x(1); ...
F/V*(u(2)-x(2))-(H/HD)*k_0*exp(-E/(R*x(2)))*x(1)-(HA/(HD*V))*(x(2)-u(3)) ...
];Объект IDNLGREY, отражающий ситуацию с моделированием, затем создается.
FileName = 'cstr_m'; % File describing the model structure. Order = [2 3 2]; % Model orders [ny nu nx]. Parameters = [1; 1; 35e6; 11850; ... % Initial parameters. 1.98589; -5960; 480; 145]; InitialStates = [8.5695; 311.267]; % Initial value of the initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, 'Name', ... 'Stirred tank reactor', ... 'TimeUnit', 'hours');
Входные параметры, состояния и выходные параметры структуры модели CSTR заданы с помощью НАБОРА методов и SETINIT. Мы также указываем, что начальные состояния по умолчанию должны быть оценены.
nlgr.InputName = {'Concentration of A in inlet feed stream' ... % u(1).
'Inlet feed stream temperature' ... % u(2).
'Jacket coolant temperature'}; % u(3).
nlgr.InputUnit = {'kgmol/m^3' 'K' 'K'};
nlgr = setinit(nlgr, 'Name', {'Concentration of A in reactor tank' ... % x(1).
'Reactor temperature'}); ... % x(2).
nlgr = setinit(nlgr, 'Unit', {'kgmol/m^3' 'K'});
nlgr = setinit(nlgr, 'Fixed', {false false});
nlgr.OutputName = {'A Concentration' ... % y(1); Concentration of A in reactor tank
'Reactor temp.'}; % y(2).
nglr.OutputUnit = {'kgmol/m^3' 'K'};
Параметры структуры модели CSTR заданы и F, V, R, и H заданы, чтобы быть зафиксированными.
nlgr = setpar(nlgr, 'Name', {'Volumetric flow rate (volume/time)' ... % F. 'Volume in reactor' ... % V. 'Pre-exponential nonthermal factor' ... % k_0. 'Activation energy' ... % E. 'Boltzmann''s ideal gas constant' ... % R. 'Heat of reaction' ... % H. 'Heat capacity times density' ... % HD. 'Overall heat transfer coefficient times tank area'}); ... % HA. nlgr = setpar(nlgr, 'Unit', {'m^3/h' 'm^3' '1/h' 'kcal/kgmol' 'kcal/(kgmol*K)' ... 'kcal/kgmol' 'kcal/(m^3*K)' 'kcal/(K*h)'}); nlgr.Parameters(1).Fixed = true; % Fix F. nlgr.Parameters(2).Fixed = true; % Fix V. nlgr.Parameters(5).Fixed = true; % Fix R. nlgr.Parameters(6).Fixed = true; % Fix H.
Посредством физического обоснования мы также знаем, что все кроме тепла параметра реакции (всегда отрицательный, потому что реакция является экзотермической) положительны. Давайте также включим это (несколько грубое) знание в нашу структуру модели CSTR:
nlgr.Parameters(1).Minimum = 0; % F. nlgr.Parameters(2).Minimum = 0; % V. nlgr.Parameters(3).Minimum = 0; % k_0. nlgr.Parameters(4).Minimum = 0; % E. nlgr.Parameters(5).Minimum = 0; % R. nlgr.Parameters(6).Maximum = 0; % H. nlgr.Parameters(7).Minimum = 0; % HD. nlgr.Parameters(8).Minimum = 0; % HA.
Сводные данные вводимой структуры модели CSTR затем получены посредством команды PRESENT:
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'cstr_m' (MATLAB file):
dx/dt = F(t, u(t), x(t), p1, ..., p8)
y(t) = H(t, u(t), x(t), p1, ..., p8) + e(t)
with 3 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 8).
Inputs:
u(1) Concentration of A in inlet feed stream(t) [kgmol/m^3]
u(2) Inlet feed stream temperature(t) [K]
u(3) Jacket coolant temperature(t) [K]
States: Initial value
x(1) Concentration of A in reactor tank(t) [kgmol/m^3] xinit@exp1 8.5695 (estimated) in [-Inf, Inf]
x(2) Reactor temperature(t) [K] xinit@exp1 311.267 (estimated) in [-Inf, Inf]
Outputs:
y(1) A Concentration(t)
y(2) Reactor temp.(t)
Parameters: Value
p1 Volumetric flow rate (volume/time) [m^3/h] 1 (fixed) in [0, Inf]
p2 Volume in reactor [m^3] 1 (fixed) in [0, Inf]
p3 Pre-exponential nonthermal factor [1/h] 3.5e+07 (estimated) in [0, Inf]
p4 Activation energy [kcal/kgmol] 11850 (estimated) in [0, Inf]
p5 Boltzmann's ideal gas constant [kcal/(kgmo..] 1.98589 (fixed) in [0, Inf]
p6 Heat of reaction [kcal/kgmol] -5960 (fixed) in [-Inf, 0]
p7 Heat capacity times density [kcal/(m^3*K)] 480 (estimated) in [0, Inf]
p8 Overall heat transfer coefficient times tank area [kcal/(K*h)] 145 (estimated) in [0, Inf]
Name: Stirred tank reactor
Status:
Created by direct construction or transformation. Not estimated.
More information in model's "Report" property.
Данные ввода - вывода
Проект эксперимента системы идентификации для многих нелинейных систем обычно намного более включается, чем для линейных систем. Это также верно для CSTR, где с одной стороны желательно, чтобы управляемый вход u_3 был таков, что это волнует систему достаточно, и с другой стороны это должно быть выбрано, чтобы быть "благоприятным для объекта" (химический процесс должен сохраниться стабильность, длительность теста должна быть максимально короткой, чтобы влиять на производство меньше всего, и т.д). Статья в [1] обсуждает выбор входных сигналов к CSTR. Там утверждается, что мультисинусоидальный вход u_3 выгоден для многоуровневого псевдослучайного входного сигнала по нескольким причинам. В идентификационных экспериментах ниже мы будем использовать два таких входных сигнала, один для оценки и один в целях валидации, которые были сгенерированы через инструмент генерации входных данных MATLAB® (графический интерфейс пользователя) любезно обеспеченный авторами упомянутой статьи.
Мы загружаем эти данные CSTR и помещаем их в два различных объекта IDDATA, ze для оценки и zv в целях валидации:
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'cstrdata')); Ts = 0.1; % 10 samples per hour! ze = iddata(y1, u1, 0.1, 'Name', 'Estimation data'); ze.InputName = nlgr.InputName; ze.InputUnit = nlgr.InputUnit; ze.OutputName = nlgr.OutputName; ze.OutputUnit = nlgr.OutputUnit; ze.Tstart = 0; ze.TimeUnit = 'hour'; ze.ExperimentName = 'Estimation'; zv = iddata(y2, u2, 0.1, 'Name', 'Validation data'); zv.InputName = nlgr.InputName; zv.InputUnit = nlgr.InputUnit; zv.OutputName = nlgr.OutputName; zv.OutputUnit = nlgr.OutputUnit; zv.Tstart = 0; zv.TimeUnit = 'hour'; zv.ExperimentName = 'Validation';
Вводы и выводы набора данных оценки ze показывают в двух графиках.
figure('Name', [ze.Name ': input data']); for i = 1:ze.Nu subplot(ze.Nu, 1, i); plot(ze.SamplingInstants, ze.InputData(:,i)); title(['Input #' num2str(i) ': ' ze.InputName{i}]); xlabel(''); axis tight; end xlabel([ze.Domain ' (' ze.TimeUnit ')']);
Рисунок 2: набор данных Оценки вводит к CSTR.
figure('Name', [ze.Name ': output data']); for i = 1:ze.Ny subplot(ze.Ny, 1, i); plot(ze.SamplingInstants, ze.OutputData(:,i)); title(['Output #' num2str(i) ': ' ze.OutputName{i}]); xlabel(''); axis tight; end xlabel([ze.Domain ' (' ze.TimeUnit ')']);
Рисунок 3: набор данных Оценки выходные параметры от CSTR.
Прежде чем мы возобновим идентификационные эксперименты, мы должны упомянуть, что сгенерированные входные сигналы обеспечивают выходные параметры CSTR, чтобы отобразить большую реакторную нелинейность (с изменениями температуры приблизительно 80 градусов в целом и того, чтобы заставлять некоторые явления "воспламенения" реактора быть очевидными). Принимая во внимание, что это волнует реактор (обычно хороший с идентификационной точки зрения), это, вероятно, а не то, как, в котором инженеры хотели бы управлять реальным реактором, особенно не тот, который является столь же экзотермическим как этот. Используя инструкции, описанные в Брауне и др. (2002), можно было затем перепроектировать эксперимент, прежде чем он будет на самом деле выполнен. В этом случае было бы интересно попытаться уменьшать длительность эксперимента и использовать мультисинусоидальные входные сигналы с более короткими длинами цикла. Цель состоит в том, чтобы уменьшать низкочастотное содержимое управляемого входного сигнала, чтобы уменьшать изменения реакторных выходных параметров.
Эффективность первоначальной модели CSTR
Насколько хороший первоначальная модель CSTR? Давайте исследуем это путем симуляции первоначальной модели с помощью входных сигналов ze и zv и давайте сравним вычисленные выходные параметры с истинными выходными параметрами (полученный путем симуляции вышеупомянутой модели IDNLGREY с помощью других параметров и добавляя некоторый шум) содержавшийся в ze и zv, соответственно.
clf compare(ze, nlgr); figure; compare(zv, nlgr);
Рисунок 4: Сравнение между истинными выходными параметрами и симулированными выходными параметрами первоначальной модели CSTR.
Оценка параметра
Соглашение между истинными и симулированными выходными параметрами первоначальной модели CSTR достойно. Чтобы далее улучшить его, мы оцениваем четыре свободных параметра модели, а также вектор начального состояния из модели при помощи набора данных оценки ze. Мы даем NLGREYEST команду отображать информацию об итерации и выполнять самое большее 25 поисковых итераций.
opt = nlgreyestOptions('Display', 'on'); opt.SearchOptions.MaxIterations = 25; nlgr = nlgreyest(ze, nlgr, opt);
При желании поиск может всегда продолжаться через второй вызов NLGREYEST (или PEM). На этот раз NLGREYEST проинструктирован, чтобы не отобразить любую информацию об итерации и только выполнить самое большее еще 5 итераций.
opt.Display = 'off';
opt.SearchOptions.MaxIterations = 5;
nlgr = nlgreyest(ze, nlgr, opt);
Эффективность предполагаемой модели CSTR
Чтобы оценить эффективность предполагаемой модели, мы еще раз используем COMPARE:
compare(ze, nlgr); figure; compare(zv, nlgr);
Рисунок 5: Сравнение между истинными выходными параметрами и симулированными выходными параметрами предполагаемой модели CSTR.
Визуальный осмотр сразу показывает, что выходные параметры предполагаемой модели близко к истинным выходным параметрам, и для ze и для zv. Улучшение является особенно значительным для набора данных валидации, где подгонки модели увеличились от отрицательных величин приблизительно до 70% и 99%, соответственно, для двух выходных параметров модели.
Дополнительная информация о предполагаемой модели CSTR затем возвращена PRESENT:
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'cstr_m' (MATLAB file):
dx/dt = F(t, u(t), x(t), p1, ..., p8)
y(t) = H(t, u(t), x(t), p1, ..., p8) + e(t)
with 3 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 8).
Inputs:
u(1) Concentration of A in inlet feed stream(t) [kgmol/m^3]
u(2) Inlet feed stream temperature(t) [K]
u(3) Jacket coolant temperature(t) [K]
States: Initial value
x(1) Concentration of A in reactor tank(t) [kgmol/m^3] xinit@exp1 8.62914 (estimated) in [-Inf, Inf]
x(2) Reactor temperature(t) [K] xinit@exp1 311.215 (estimated) in [-Inf, Inf]
Outputs:
y(1) A Concentration(t)
y(2) Reactor temp.(t)
Parameters: Value Standard Deviation
p1 Volumetric flow rate (volume/time) [m^3/h] 1 0 (fixed) in [0, Inf]
p2 Volume in reactor [m^3] 1 0 (fixed) in [0, Inf]
p3 Pre-exponential nonthermal factor [1/h] 3.55889e+07 17548.3 (estimated) in [0, Inf]
p4 Activation energy [kcal/kgmol] 11853.9 0.0703052 (estimated) in [0, Inf]
p5 Boltzmann's ideal gas constant [kcal/(kgmo..] 1.98589 0 (fixed) in [0, Inf]
p6 Heat of reaction [kcal/kgmol] -5960 0 (fixed) in [-Inf, 0]
p7 Heat capacity times density [kcal/(m^3*K)] 500.71 0.194139 (estimated) in [0, Inf]
p8 Overall heat transfer coefficient times tank area [kcal/(K*h)] 150.127 0.0697455 (estimated) in [0, Inf]
Name: Stirred tank reactor
Status:
Termination condition: Maximum number of iterations or number of function evaluations reached..
Number of iterations: 6, Number of function evaluations: 7
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Estimation data".
Fit to estimation data: [71.36;99.18]%
FPE: 0.02736, MSE: 0.6684
More information in model's "Report" property.
Симуляция предполагаемой модели CSTR в Simulink
Модель IDNLGREY может также импортироваться и использоваться в Simulink. Модель Simulink "cstr_sim" импортирует набор данных валидации zv и передает его данные блоку модели IDNLGREY Simulink, который, когда симулировано производит выходные параметры, которые вместе с используемыми входными сигналами хранятся в рабочем пространстве MATLAB в объекте IDDATA zsim. (Последние пять строк кода ниже используются, чтобы гарантировать, что соответствующие входные параметры модели питаются модель Simulink. Это только необходимо, если idnlgreydemo9 запущен в функции, где доступ к zv и nlgr не может быть гарантирован.)
open_system('cstr_sim'); if ~evalin('base', 'exist(''zv'', ''var'')') cstrws = get_param(bdroot, 'modelworkspace'); cstrws.assignin('zv', zv); cstrws.assignin('nlgr', nlgr); end
Рисунок 6: модель Simulink, содержащая предполагаемую модель CSTR.
Типовой блок Библиотеки Simulink IDNLGREY найден в стандартной Библиотеке Simulink системы идентификации и может копироваться в и использоваться в любой модели Simulink. Например, в случае CSTR это могло очень хорошо использоваться в расположении управления с обратной связью.
Блок IDNLGREY должен быть сконфигурирован, прежде чем он будет симулирован. Это сделано путем ввода переменной рабочего пространства MATLAB содержание модели IDNLGREY (здесь nlgr), или путем определения соответствующего объекта модели IDNLGREY с помощью idnlgrey конструктора в окне редактирования параметра пометил "IDNLGREY model". Здесь также возможно задать вектор начального состояния, чтобы использовать (значением по умолчанию является внутренне сохраненный вектор начального состояния из заданного объекта IDNLGREY).
Объект модели IDNLGREY хранит свойства, задающие настройку решателя дифференциала/разностного уравнения, используемого СИМОМ, PREDICT, NLGREYEST и так далее (см. nlgr.SimulationOptions для опций, управляющих симуляцией модели). В Simulink всегда заменяются эти настройки так, чтобы опции Simulink указали, что решатель используется. Если объект модели IDNLGREY задает и использует другой решатель, то результат симуляции, полученный в Simulink, может отличаться от полученного с IDNLGREY/SIM в MATLAB. Пример, иллюстрирующий это, обеспечивается в примере "idnlgreydemo10".
С улаженными опциями решателя мы можем затем выполнить симуляцию командной строки cstr_sim модели Simulink (который здесь будет проводиться для предполагаемой модели CSTR). (Вызов evalin необходим, чтобы получить zsim в случае, если idnlgreydemo9 запущен в функции.)
sim('cstr_sim'); if ~exist('zsim', 'var') zsim = evalin('base', 'zsim'); end zsim.InputName = nlgr.InputName; zsim.InputUnit = nlgr.InputUnit; zsim.OutputName = nlgr.OutputName; zsim.OutputUnit = nlgr.OutputUnit; zsim.TimeUnit = 'hour';
Давайте наконец придем к заключению, что пример путем графического вывода Simulink получил результат симуляции.
figure('Name', 'Simulink simulation result of estimated CSTR model'); for i = 1:zsim.Ny subplot(zsim.Ny, 1, i); plot(zsim.SamplingInstants, ze.OutputData(:,i)); title(['Output #' num2str(i) ': ' zsim.OutputName{i}]); xlabel(''); axis tight; end xlabel([zsim.Domain ' (' zsim.TimeUnit ')']);
Рисунок 7: Выходные параметры, полученные путем симуляции предполагаемой модели CSTR в Simulink.
Заключения
Этот пример касался моделирования и идентификации неадиабатического непрерывного реактора смесителя. В частности, это было проиллюстрировано, как импортировать и использовать модель IDNLGREY в Simulink.
Ссылки
[1] Браун, M.W., Р. Ортис-Мохика и Д.Е. Ривера, "Приложение минимального фактора гребня мультисинусоидальные сигналы для 'благоприятной для объекта' идентификации нелинейных систем процесса", в Практике Разработки Управления, Издании 10, № 3, 2002, стр 301-313.