В этом примере показано, как включить и смоделировать модель 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 моделируется с использованием основных принципов учета и энергосбережения. Изменение концентрации реагента А в сосуде в единицу времени 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 член теплоемкости, член плотности, 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 и энергию Е активации в экспериментах со стендовой шкалой. Это упростило бы задачу идентификации, чтобы рассмотреть только два неизвестных: теплоемкость 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)Эта информация вводится в файл 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? Рассмотрим это путем моделирования исходной модели с использованием входных сигналов 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);
Чтобы оценить производительность расчетной модели, мы снова используем COMPARE:
compare(ze, nlgr); figure; compare(zv, nlgr);
Рис. 5: Сравнение между истинными выходами и моделируемыми выходами расчетной модели CSTR.
Визуальный осмотр сразу показывает, что выходы расчетной модели близки к истинным выходам, как для ze, так и для zv. Улучшение особенно важно для набора проверочных данных, где посадки модели увеличились с отрицательных значений до около 70% и 99% соответственно для двух выходов модели.
Следующая информация о предполагаемой модели CSTR возвращается в разделе НАСТОЯЩЕЕ:
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.
Модель 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 model». Здесь также можно указать используемый вектор начального состояния (по умолчанию - внутренний сохраненный вектор начального состояния указанного объекта IDNLGREY).
Объект модели IDNLGREY хранит свойства, определяющие настройку решателя дифференциальных/разностных уравнений, используемого SIM, PREDICT, NLGREYEST и т.д. (см. nlgr. Опции для опций, управляющих моделированием модели). В 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] Браун, М.В., Р. Ортис-Мохика и Д.Е. Ривера, «Применение мультисинусоидальных сигналов минимального коэффициента гребня для» удобной для установки «идентификации нелинейных технологических систем», в Практике проектирования управления, том 10, № 3, 2002, стр. 301-313.