В этом примере показано, как выполнить IDNLGREY, моделирующий на основе файлов модели C MEX. Это использует простую систему, где нелинейное пространство состояний, моделирующее действительно, окупается.
Цель состоит в том, чтобы смоделировать уровень жидкости более низкого бака лабораторных весов две системы бака, схематично показанные в рисунке 1.
Рисунок 1: Схематическое представление двух систем бака.
Мы запускаем задание моделирования путем загрузки доступных данных ввода - вывода, которые были симулированы с помощью ниже структуры модели IDNLGREY с шумом, добавленным к выходу. twotankdata.mat файл содержит один набор данных с 3 000 выборок ввода - вывода, сгенерированного использования уровня выборки (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)
где Ай [м^2] является площадью поперечного сечения бака i и Qini (t), и Qouti (t) [м^3/c] являются притоком в и оттоком от бака 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 на C 6 параметрами (или константы), A1, k, a1, g, A2 и a2. Файл MEX на C обычно немного более включается, чем соответствующий файл записанный язык MATLAB использования, но MEX C, моделирующий обычно, дает явное преимущество в терминах скорости выполнения, специально для более сложных моделей. Шаблон файл MEX на C обеспечивается (см. ниже) помочь пользователю структурировать код. Для большинства приложений это достаточно, чтобы задать количество выходных параметров и ввести строки кода, которые описывают дуплекс и 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 скалярных параметров, и следовательно количество входных параметров к файлу моделирования MEX C должно быть 3+Npo = 3+6 = 9. Запаздывание 10:th аргумент может здесь быть не использован, когда никакой дополнительный FileArgument не нанимается в этом приложении.
Записи файла моделирования MEX C обычно делаются на четырех шагах:
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: C исходный код MEX для двух систем бака.
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. Обе этих функции содержат списки аргументов с выходом, который будет вычислен (дуплекс или y) в положении 1, после которого следует за всеми переменными и параметрами, требуемыми вычислить правую сторону (стороны) состояния и выходных уравнений, соответственно.
Первый шаг в этих функциях должен распаковать параметры модели, которые будут использоваться в последующих уравнениях. Любое допустимое имя переменной (за исключением используемых в списке входных параметров) может использоваться, чтобы обеспечить физически понятные имена отдельных параметров.
Как имеет место в C, первый элемент массива хранится в положении 0. Следовательно, дуплекс [0] в C соответствует дуплексу (1) в MATLAB® (или только дуплекс в случае, если это - скаляр), вход u [0] соответствует вам (или u (1)), параметр A1 [0] соответствует A1 и так далее.
Два файла модели бака включают расчеты квадратного корня. Это включено посредством включения математической библиотеки C math.h. Математическая библиотека осознает наиболее распространенные тригонометрические функции (sin, потому что, tan, asin, acos, atan, и т.д.), экспоненциал (exp) и логарифмы (журнал, log10), квадратный корень (sqrt) и степень функций (голова) и расчеты абсолютного значения (fabs). math.h библиотека должна быть включена каждый раз, когда любая функция math.h используется; в противном случае это может быть не использовано. См. "Примеры на Нелинейной Серой Идентификации Модели Поля: Создание Файлы Модели IDNLGREY" для получения дальнейшей информации о математической библиотеке C.
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, описывающий две системы бака. Для удобства мы также устанавливаем некоторую бухгалтерскую информацию о вводах и выводах (имя и модули).
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 мы симулируем систему с помощью начальных значений параметров. Мы используем решатель для дифференциальных уравнений по умолчанию (45 решателей Рунге-Кутта с адаптивной корректировкой длины шага) и устанавливаем погрешности абсолютной и относительной погрешности довольно маленьким значениям (1e-6 и 1e-5, соответственно). Заметьте, что команда COMPARE, когда названо двумя входными параметрами, когда значение по умолчанию оценит все начальное состояние (состояния) независимо от того, было ли какое-либо начальное состояние задано, чтобы быть 'Зафиксированным'. Для того, чтобы только оценить свободное начальное состояние (состояния), вызовите COMPARE с одной третью и четвертым входным параметром можно следующим образом: сравните (z, nlgr, 'init', 'm'); когда оба начальных состояния модели бака по умолчанию 'Фиксируются', никакая оценка начального состояния не будет выполняться этой командой.
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 В пошаговым способом. Мы делаем это путем вызова ШАГА с различными заданными амплитудами шага, начинающими с фиксированного смещения 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.