Система с двумя баками: C Файл MEX на C Моделирование непрерывной по времени системы SISO
В этом примере показано, как выполнить моделирование IDNLGREY на основе файлов модели C MEX. Он использует простую систему, где нелинейное моделирование пространства состояний действительно окупается.
Система с двумя баками
Цель состоит в том, чтобы смоделировать уровень жидкости в нижнем баке лаборатории шкалы двух баков систем, схематически показанных на Фигуру 1.
Фигура 1: схематический вид двух баков систем.
Входные-выходные данные
Мы начинаем задание моделирования с загрузки доступных входно-выходных данных, которые были симулированы с использованием структуры модели IDNLGREY ниже, с добавлением шума к выходу. Файл twotankdata.mat содержит один набор данных с 3000 входно-выходных выборок, сгенерированных с использованием частоты дискретизации (Ts) 0,2 секунды. Вход u (t) является напряжением [V], приложенным к насосу, который генерирует приток к верхнему баку. Довольно небольшое отверстие в нижней части этого верхнего бака приводит к выходу, который поступает в нижний бак, и выход y (t) системы двух баков затем является уровнем жидкости [ m ] нижнего бака. Создадим объект IDDATA z, чтобы сохранить данные бака. В целях бухгалтерского учета и документации мы также указываем имена каналов и модулей. Этот шаг является необязательным.
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'twotankdata')); z = iddata(y, u, 0.2, 'Name', 'Two tanks'); set(z, 'InputName', 'Pump voltage', 'InputUnit', 'V', ... 'OutputName', 'Lower tank water level', 'OutputUnit', 'm', ... 'Tstart', 0, 'TimeUnit', 's');
Входно-выходные данные, которые будут использоваться для оценки, показаны в окне графика.
figure('Name', [z.Name ': input-output data']); plot(z);
Фигура 2: Входно-выходные данные из системы с двумя баками.
Моделирование системы с двумя баками
Следующим шагом является определение структуры модели, описывающей две системы баков. Для этого пусть x1 (t) и x2 (t) обозначают уровень воды в верхнем и нижнем баке. Для каждого бака фундаментальная физика (баланс масс) утверждает, что изменение объема воды зависит от различия между входным и выходным потоками как (i = 1, 2):
d/dt (Ai*xi(t)) = Qini(t) - Qouti(t)
где Ai [m ^ 2] - площадь поперечного сечения бака i, а Qini (t) и Qouti (t) [m ^ 3/s] - приток и отток из бака i в момент t.
Для верхнего бака приток принимается пропорциональным напряжению, приложенному к насосу, т.е. Qin1 (t) = k * u (t). Поскольку выходное отверстие верхнего бака мало, может быть применен закон Бернулли, утверждающий, что поток пропорционален квадратному корню уровня воды или, точнее:
Qout1(t) = a1*sqrt(2*g*x1(t))
где a1 - площадь поперечного сечения выходного отверстия, а g - гравитационная константа. Для нижнего бака приток равен оттоку из верхнего бака, т.е. Qin2 (t) = Qout1 (t), и отток дается законом Бернулли:
Qout2(t) = a2*sqrt(2*g*x2(t))
где a2 - площадь поперечного сечения выходного отверстия.
В целом эти факты приводят к следующей структуре в пространстве состояний:
d/dt x1(t) = 1/A1*(k*u(t) - a1*sqrt(2*g*x1(t)))
d/dt x2(t) = 1/A2*(a1*sqrt(2*g*x1(t)) - a2*sqrt(2*g*x2(t)))
y(t) = x2(t)Два файла модели MEX Tank C
Эти уравнения затем помещаются в файл MEX на C с 6 параметрами (или константами), A1, k, a1, g, A2 и a2. Обычно файл MEX на C является битом более вовлеченным, чем соответствующий файл, написанный с использованием языка MATLAB, но моделирование C MEX обычно дает отчетливое преимущество с точки зрения скорости выполнения, особенно для более сложных моделей. Для помощи пользователю в структурировании кода предоставляется шаблон С файла MEX на C (см. ниже). Для большинства приложений достаточно задать количество выходов и ввести в этот шаблон кода линий, которые описывают dx и y. Файл MEX на C IDNLGREY всегда должна быть структурирована, чтобы вернуть два выхода:
dx: the right-hand side(s) of the state-space equation(s) y: the right-hand side(s) of the output equation(s)
и это должно занять 3 + Npo (+ 1) входных параметров, заданных следующим образом:
t: the current time
x: the state vector at time t ([] for static models)
u: the input vector at time t ([] for time-series models)
p1, p2, ..., pNpo: the individual parameters (which can be real
scalars, column vectors or 2-dimensional matrices); Npo is here
the number of parameter objects, which for models with scalar
parameters coincide with the number of parameters Np
FileArgument: optional inputs to the model fileВ нашей двух баков системе существует 6 скалярных параметров, и, следовательно, количество входных параметров в файл моделирования C MEX должно быть 3 + Npo = 3 + 6 = 9. Конечный аргумент 10: th здесь может быть опущен, так как в этом приложении не используется необязательный FileArgument.
Запись файла моделирования C MEX обычно выполняется в четыре этапа:
1. Inclusion of C-libraries and definitions of the number of outputs.
2. Writing the function computing the right-hand side(s) of the state
equation(s), compute_dx.
3. Writing the function computing the right-hand side(s) of the output
equation(s), compute_y.
4. Writing the main interface function, which includes basic error
checking functionality, code for creating and handling input and
output arguments, and calls to compute_dx and compute_y.Let us view the C MEX source file (except for some comments) for the two tank system and based on this discuss these four items in some more detail.
Фигура. 3. Исходный код MEX на C для двух резервуарных систем.
1. Две библиотеки C mex.h и math.h обычно включены, чтобы предоставить доступ к ряду связанных с MEX функций, а также математических функций. Количество выходов также объявляется на каждый файл моделирования с помощью стандартного C-define:
/* Include libraries. */ #include "mex.h" #include "math.h"
/* Specify the number of outputs here. */ #define NY 1
2-3. Далее в файле мы находим функции для обновления состояний, compute_dx и выхода, compute_y. Обе эти функции содержат списки аргументов с выходом, который будет вычислен (dx или y) в позиции 1, после чего следуют всем переменным и параметрам, необходимым для вычисления правой стороны (ей) состояния и выходных уравнений, соответственно.
Первый шаг в этих функциях - распаковать параметры модели, которые будут использоваться в последующих уравнениях. Любое допустимое имя переменной (кроме тех, что используются в списке входных параметров) может использоваться, чтобы предоставить физически значимые имена отдельных параметров.
Как и в случае С, первый элемент массива сохранен в положении 0. Следовательно, dx [0] в C соответствует dx (1) в MATLAB ® (или просто dx в случае, если это скаляр), вход u [0] соответствует u (или u (1)), параметр A1 [0] соответствует A1 и так далее.
Два файла модели бака включают квадратные корневые расчеты. Это включено путем включения математической библиотеки C math.h. Математическая библиотека реализует наиболее распространенные тригонометрические функции (sin, cos, tan, asin, acos, atan и т. д.), экспоненциал (exp) и логарифмы (журнал, log10), квадратный корень (sqrt) и степень функций (pow), и абсолютные расчеты значений (fabs). Библиотека math.h должна включаться всякий раз, когда используется какая-либо функция math.h; в противном случае его можно исключить. Для получения дополнительной информации о математической библиотеке C смотрите «Руководства по идентификации модели нелинейного серого ящика: Создание файлов модели IDNLGREY».
4. Основная функция интерфейса должна почти всегда иметь одинаковое содержимое и для большинства приложений никакая модификация не нужна. В принципе, единственная часть, которая может быть рассмотрена для изменений, это место, где выполняются вызовы compute_dx и compute_y. Для статических систем можно исключить вызов, чтобы compute_dx. В других ситуациях может быть желательно только передать переменные и параметры, указанные в уравнениях состояния и выходе. Например, в выходном уравнении двух баков, где используется только одно состояние, можно было бы очень сократить список входных параметров до:
void compute_y(double *y, double *x)
и вызывайте compute_y как:
compute_y(y, x);
Входные параметры compute_dx и compute_y могут также быть расширены, чтобы включать дополнительные переменные, выведенные в функции интерфейса, такие как количество состояний и количество параметров.
После того, как исходный файл модели будет завершен, он должен быть скомпилирован, что можно сделать из командной строки MATLAB с помощью команды mex; см. «Помощь mex». (Этот шаг здесь опущен.)
При разработке конкретных моделей Файлов MEX на C часто полезно начать работу, скопировав файл шаблона IDNLGREY C MEX. Этот шаблон содержит исходный код скелета, а также подробные инструкции по настройке кода для конкретного приложения. Расположение файла шаблона отображается путем ввода следующего в командной строке MATLAB.
fullfile(matlabroot, 'toolbox', 'ident', 'nlident', 'IDNLGREY_MODEL_TEMPLATE.c')
Также смотрите пример «Создание файлов модели IDNLGREY» для получения дополнительной информации о файлах модели IDNLGREY C MEX.
Создание объекта модели IDNLGREY двух баков
Следующим шагом является создание объекта IDNLGREY, описывающего систему двух баков. Для удобства мы также устанавливаем некоторую бухгалтерскую информацию о входах и выходах (имя и модули).
FileName = 'twotanks_c'; % File describing the model structure. Order = [1 1 2]; % Model orders [ny nu nx]. Parameters = {0.5; 0.0035; 0.019; ... 9.81; 0.25; 0.016}; % Initial parameters. InitialStates = [0; 0.1]; % Initial value of initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Two tanks'); set(nlgr, 'InputName', 'Pump voltage', 'InputUnit', 'V', ... 'OutputName', 'Lower tank water level', 'OutputUnit', 'm', ... ... 'TimeUnit', 's');
Мы продолжаем добавлять информацию об именах и модулях состояний и параметрах модели через команды SETINIT и SETPAR. Кроме того, оба состояния x1 (t) и x2 (t) являются уровнями бака, которые не могут быть отрицательными, и, таким образом, мы также задаем, что x1 (0) и x2 (0) > = 0 через свойство 'Minimum'. На самом деле, мы также знаем, что все параметры модели должны быть строго положительными. Поэтому мы устанавливаем свойство 'Minimum' всех параметров на некоторое маленькое положительное значение (eps (0)). Эти настройки подразумевают, что оценка ограничений будет выполнена на предстоящем шаге оценки (т.е. предполагаемая модель будет такой моделью, чтобы все введенные ограничения выполнялись).
nlgr = setinit(nlgr, 'Name', {'Upper tank water level' 'Lower tank water level'}); nlgr = setinit(nlgr, 'Unit', {'m' 'm'}); nlgr = setinit(nlgr, 'Minimum', {0 0}); % Positive levels! nlgr = setpar(nlgr, 'Name', {'Upper tank area' ... 'Pump constant' ... 'Upper tank outlet area' ... 'Gravity constant' ... 'Lower tank area' ... 'Lower tank outlet area'}); nlgr = setpar(nlgr, 'Unit', {'m^2' 'm^3/(s*V)' 'm^2' 'm/(s^2)' 'm^2' 'm^2'}); nlgr = setpar(nlgr, 'Minimum', num2cell(eps(0)*ones(6,1))); % All parameters > 0!
Площади поперечных сечений ( A1 и A2) двух баков можно определить довольно точно. Поэтому мы рассматриваем эти и g как константы и проверяем, что поле 'Fixed' правильно установлено для всех 6 параметров через команду GETPAR. В целом это означает, что будут оценены 3 из параметров модели.
nlgr.Parameters(1).Fixed = true;
nlgr.Parameters(4).Fixed = true;
nlgr.Parameters(5).Fixed = true;
getpar(nlgr, 'Fixed')
ans =
6x1 cell array
{[1]}
{[0]}
{[0]}
{[1]}
{[1]}
{[0]}
Эффективность начальной модели двух баков
Перед оценкой свободных параметров k, a1 и a2 моделируем систему с помощью начальных значений параметров. Мы используем решатель для дифференциальных уравнений по умолчанию (решатель Runge-Kutta 45 с адаптивной регулировкой длины шага) и устанавливаем абсолютные и относительные допуски на довольно маленькие значения (1e-6 и 1e-5, соответственно). Заметьте, что команда COMPARE, при вызове с двумя входными параметрами, по умолчанию будет оценивать все начальные состояния (состояния ) независимо от того, определено ли какое-либо начальное состояние как 'Fixed'. В порядок только для оценки свободных начальных состояний (состояний ), вызова COMPARE с третьим и четвертым входным параметром следующим образом: compare (z, nlgr, 'init', 'm'); поскольку оба начальных состояния модели бака по умолчанию являются 'Fixed', никакая начальная оценка состояния не будет выполнена этой командой.
nlgr.SimulationOptions.AbsTol = 1e-6; nlgr.SimulationOptions.RelTol = 1e-5; compare(z, nlgr);
Фигура 4: Сравнение истинного выхода и моделируемого выхода исходной модели двух баков.
Моделируемые и истинные выходы показаны в окне графика, и, как видно, подгонка не так впечатляет.
Оценка параметра
В порядок улучшения подгонки следующие 3 свободных параметра оцениваются с помощью NLGREYEST. (Поскольку по умолчанию поля 'Fixed' всех начальных состояний являются ложными, оценка начальных состояний не будет выполнена в этом вызове оценщику.)
nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on'));
Эффективность предполагаемой модели двух баков
Чтобы исследовать эффективность предполагаемой модели, выполняется её симуляция (начальные состояния здесь переоцениваются).
compare(z, nlgr);
Фигура 5: Сравнение истинного выхода и моделируемого выхода предполагаемой модели двух баков.
Соглашение между истинным и моделируемым выходами довольно хорошее. Остающийся вопрос, однако, заключается в том, могут ли две системы баков быть точно описаны с помощью более простой и линейной структуры модели. Чтобы ответить на это, давайте попробуем подогнать данные к некоторым стандартным линейным структурам модели, а затем используйте COMPARE, чтобы увидеть, насколько хорошо эти модели захватывают динамику баков.
nk = delayest(z); arx22 = arx(z, [2 2 nk]); % Second order linear ARX model. arx33 = arx(z, [3 3 nk]); % Third order linear ARX model. arx44 = arx(z, [4 4 nk]); % Fourth order linear ARX model. oe22 = oe(z, [2 2 nk]); % Second order linear OE model. oe33 = oe(z, [3 3 nk]); % Third order linear OE model. oe44 = oe(z, [4 4 nk]); % Fourth order linear OE model. sslin = ssest(z); % State-space model (order determined automatically) compare(z, nlgr, 'b', arx22, 'm-', arx33, 'm:', arx44, 'm--', ... oe22, 'g-', oe33, 'g:', oe44, 'g--', sslin, 'r-');
Фигура. 6. Сравнение истинного выхода и моделируемых выходов ряда предполагаемых двух моделей бака.
График сравнения ясно показывает, что линейные модели не могут подобрать всю динамику двух баков. Предполагаемая нелинейная модель IDNLGREY, с другой стороны, показывает превосходную подгонку истинному выходу. В сложение параметры модели IDNLGREY также хорошо линия с параметрами, используемыми для генерации истинного выхода. В следующих расчетах отображения мы используем команду GETPVEC, которая возвращает вектор параметра, созданный из массива структур, содержащего параметры модели объекта IDNLGREY.
disp(' True Estimated parameter vector'); ptrue = [0.5; 0.005; 0.02; 9.81; 0.25; 0.015]; fprintf(' %1.4f %1.4f\n', [ptrue'; getpvec(nlgr)']);
True Estimated parameter vector 0.5000 0.5000 0.0050 0.0049 0.0200 0.0200 9.8100 9.8100 0.2500 0.2500 0.0150 0.0147
Ошибки предсказания, полученные с помощью PE, малы и очень напоминают случайный шум.
figure; pe(z, nlgr);
Фигура 7: Ошибки предсказания, полученные для предполагаемой модели IDNLGREY с двумя баками.
Давайте также исследуем, что произойдет, если входное напряжение будет увеличено с 5 до 6, 7, 8, 9 и 10 В пошаговым образом. Мы делаем это, вызывая STEP с различными заданными амплитудами шага, начиная с фиксированного смещения 5 В. Переходная характеристика строения облегчается выделенным набором параметров, созданным stepDataOptions:
figure('Name', [nlgr.Name ': step responses']); t = (-20:0.1:80)'; Opt = stepDataOptions('InputOffset',5,'StepAmplitude',6); step(nlgr, t, 'b', Opt); hold on Opt.StepAmplitude = 7; step(nlgr, t, 'g', Opt); Opt.StepAmplitude = 8; step(nlgr, t, 'r', Opt); Opt.StepAmplitude = 9; step(nlgr, t, 'm', Opt); Opt.StepAmplitude = 10; step(nlgr, t, 'k', Opt); grid on; legend('5 -> 6 V', '5 -> 7 V', '5 -> 8 V', '5 -> 9 V', '5 -> 10 V', ... 'Location', 'Best');
Фигура 8: Переходные характеристики, полученные для расчетной модели IDNLGREY с двумя баками.
Наконец, используя команду PRESENT, мы получаем сводную информацию о предполагаемой модели:
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'twotanks_c' (MEX-file):
dx/dt = F(t, u(t), x(t), p1, ..., p6)
y(t) = H(t, u(t), x(t), p1, ..., p6) + e(t)
with 1 input(s), 2 state(s), 1 output(s), and 3 free parameter(s) (out of 6).
Inputs:
u(1) Pump voltage(t) [V]
States: Initial value
x(1) Upper tank water level(t) [m] xinit@exp1 0 (fixed) in [0, Inf]
x(2) Lower tank water level(t) [m] xinit@exp1 0.1 (fixed) in [0, Inf]
Outputs:
y(1) Lower tank water level(t) [m]
Parameters: Value Standard Deviation
p1 Upper tank area [m^2] 0.5 0 (fixed) in ]0, Inf]
p2 Pump constant [m^3/(s*V)] 0.00488584 0.0259032 (estimated) in ]0, Inf]
p3 Upper tank outlet area [m^2] 0.0199719 0.0064682 (estimated) in ]0, Inf]
p4 Gravity constant [m/(s^2)] 9.81 0 (fixed) in ]0, Inf]
p5 Lower tank area [m^2] 0.25 0 (fixed) in ]0, Inf]
p6 Lower tank outlet area [m^2] 0.0146546 0.0776058 (estimated) in ]0, Inf]
Name: Two tanks
Status:
Termination condition: Change in cost was less than the specified tolerance..
Number of iterations: 8, Number of function evaluations: 9
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Two tanks".
Fit to estimation data: 97.35%
FPE: 2.419e-05, MSE: 2.414e-05
More information in model's "Report" property.
Заключения
В этом примере мы показали:
1. how to use C MEX-files for IDNLGREY modeling, and
2. provided a rather simple example where nonlinear state-space
modeling shows good potentialДополнительная информация
Для получения дополнительной информации об идентификации динамических систем с помощью System Identification Toolbox™ посетите информационную страницу System Identification Toolbox.