В этом примере показано, как выполнять моделирование 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))
где а2 - площадь поперечного сечения выпускного отверстия.
В совокупности эти факты приводят к следующей государственно-космической структуре:
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)Эти уравнения затем помещаются в C MEX-файл с 6 параметрами (или константами), A1, k, a1, g, A2 и a2. C MEX-файл обычно задействован немного больше, чем соответствующий файл, написанный на языке MATLAB, но моделирование C MEX обычно дает явное преимущество с точки зрения скорости выполнения, особенно для более сложных моделей. Шаблон C MEX-файла предоставляется (см. ниже), чтобы помочь пользователю структурировать код. Для большинства приложений достаточно определить количество выходов и ввести в этот шаблон строки кода, описывающие dx и y. Файл IDNLGREY C MEX всегда должен быть структурирован для получения двух выходных данных:
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: Исходный код 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. Обе эти функции содержат списки аргументов, при этом выходные данные вычисляются (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) и логарифмы (log, 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». (Здесь этот шаг опущен.)
При разработке специфичных для модели C MEX-файлов часто полезно начинать работу с копирования файла шаблона IDNLGREY C MEX. Этот шаблон содержит исходный код скелета, а также подробные инструкции по настройке кода для конкретного приложения. Расположение файла шаблона отображается путем ввода следующего в командной строке MATLAB.
fullfile(matlabroot, 'toolbox', 'ident', 'nlident', 'IDNLGREY_MODEL_TEMPLATE.c')
Для получения дополнительной информации о файлах моделей 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'); поскольку оба начальных состояния модели резервуара по умолчанию являются фиксированными, оценка начального состояния по этой команде не выполняется.
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: Сравнение истинных выходных данных и смоделированных выходных данных расчетной модели двух резервуаров.
Соглашение между истинными и смоделированными выходами довольно хорошее. Однако остается вопрос, можно ли точно описать систему двух резервуаров, используя более простую и линейную структуру модели. Чтобы ответить на этот вопрос, попытаемся привести данные в соответствие с некоторыми стандартными линейными модельными структурами, а затем используйте команду СРАВНИТЬ, чтобы увидеть, насколько хорошо эти модели отражают динамику танков.
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Дополнительные сведения об идентификации динамических систем с помощью Toolbox™ System Identification см. на информационной странице System Identification Toolbox.