Этот пример показывает, как выполнить 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)
где Ай [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 на 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, потому что, загар, 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.