Неадиабатический непрерывный реактор с мешалкой: Моделирование файлов MATLAB с симуляциями в Simulink
®В этом примере показано, как включить и симулировать модель IDNLGREY в Simulink ®. Мы используем химическую реакционную систему в качестве базиса для моделирования. Первая часть примера, посвященная моделированию и идентификации, может быть запущена без Simulink.
Моделирование неадиабатического непрерывного реактора с Баком перемешиванием
Довольно распространенной химической системой, с которой сталкиваются в технологической промышленности, является реактор непрерывно перемешиваемого Бака (CSTR). Здесь мы изучим сахарный диабатический (то есть неадиабатический) бак, широко описанный в книге Бекетта «Динамика процесса: моделирование, анализ и симуляция», опубликованной Prentice Hall, 1998. Сосуд принимается идеально смешанным, и происходит одна экзотермическая и необратимая реакция первого порядка A -- > B. Принципиальная схема корпуса и окружающей охлаждающей рубашки показана в окне графика. Заметьте, что это эскиз; в действительности поток хладагента обычно окружает всю рубашку реактора, а не только ее дно.
Фигура 1: Принципиальная схема CSTR.
Модель CSTR требуется для более совершенных подходов к управлению. Входной поток реагента А подают с постоянной скоростью F. После перемешивания конечный продукт поступает из емкости с той же скоростью, что и реагент А, в бак (объем 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 - энергия активации, идеальная газовая константа R Больцмана и 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 термин теплоемкости, rho a термин плотности, U полный коэффициент теплопередачи и A площадь для теплообмена (площадь теплоносителя/сосуда).
В совокупности 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].
которые также являются естественными состояниями модели, т.е. 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.Четыре из параметров здесь заданы как свободные (т.е. будут оценены). На практике, вероятно, можно было бы также определить предварительный экспоненциал нетермический фактор k_0 и энергию активации E от стендовых шкал экспериментов. Это упростит задачу идентификации только для рассмотрения двух неизвестных: теплоемкости c_p и общего коэффициента теплопередачи U (которые выводятся из HD и HA, соответственно, как известны rho и 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 задаются с помощью методов SET и 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 ® (GUI), любезно предоставленный авторами упомянутой статьи.
Мы загружаем эти данные 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 степеней в целом и вызывая некоторые явления «воспламенения» реактора, быть заметными). В то время как это возбуждает реактор (обычно хорошо с точки зрения идентификации), это, вероятно, не тот способ, которым инженеры хотели бы управлять реальным реактором, особенно не такой экзотермический, как этот. Используя рекомендации, описанные в Braun et al. (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 и передает его данные в блок модели Simulink IDNLGREY, который при моделировании производит выходы, которые вместе с используемыми входными сигналами хранятся в рабочем пространстве 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.
Типовой блок библиотеки IDNLGREY Simulink находится в стандартной библиотеке системы идентификации Simulink и может быть скопирован и использован в любой модели Simulink. Для примера, в случае CSTR его вполне можно использовать в системе управления с обратной связью.
Перед моделированием блок IDNLGREY должен быть сконфигурирован. Это достигается путем ввода переменной рабочего пространства MATLAB, содержащей модель IDNLGREY (здесь nlgr), или путем определения соответствующего объекта модели IDNLGREY с помощью конструктора idnlgrey в поле редактирования параметра, обозначенном как «модель IDNLGREY Y». Здесь также можно задать вектор начального состояния, который будет использоваться (по умолчанию это внутренне сохраненный вектор начального состояния указанного объекта IDNLGREY).
Объект модели IDNLGREY хранит свойства, задающие настройку решателя дифференциального/разностного уравнения, используемого SIM, PREDICT, NLGREYEST и так далее (см. nlgr. SimulationOptions для опций, управляющих симуляцией модели). В Simulink эти настройки всегда переопределяются, так что используются опции Simulink заданного решателя. Если объект модели IDNLGREY задает и использует другой решатель, то результат симуляции, полученный в Simulink, может отличаться от результата, полученного с помощью IDNLGREY/SIM в MATLAB. Пример, иллюстрирующий это, приведен в примере «idnlgreydemo10».
С установленными опциями решателя мы можем далее выполнить симуляцию командной строки cstr_sim модели Simulink (которая здесь будет проведена для предполагаемой модели CSTR). (Вызов оценки необходим для извлечения 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] Braun, M.W., R. Ortiz-Mojica and D.E. Rivera, «Application of minimum crest factor multisinusoidal signals for 'plant-friendly' identification of nonlinear process systems», Control.