В этом примере показано, как записать файлы ODE для нелинейных серых моделей в виде файлов MATLAB и C MEX.
Моделирование серого ящика концептуально отличается от моделирования черного ящика тем, что оно включает в себя более комплексный шаг моделирования. Для IDNLGREY (нелинейный объект модели «серый ящик»; нелинейный аналог IDGREY), этот шаг состоит из создания файла ОДУ, также называемого «файлом модели». Файл ОДУ определяет правые стороны состояния и выходные уравнения, к которым обычно приходят с помощью физического моделирования первого принципа. В этом примере мы сосредоточимся на общих аспектах его реализации в виде файла MATLAB или файла C MEX.
IDNLGREY поддерживает оценку параметров и начальных состояний в нелинейных структурах модели, написанных на следующей явной форме «состояние-пространство» (так называемая форма «вывод-ошибка», «OE», названная так, что шум e (t) влияет только на выход структуры модели аддитивным образом):
xn(t) = F(t, x(t), u(t), p1, ..., pNpo); x(0) = X0;
y(t) = H(t, x(t), u(t), p1, ..., pNpo) + e(t)
Для дискретных временных структур xn (t) = x (T + Ts), где Ts является временем выборки, и для непрерывных временных структур xn (t) = d/dt x (t). Кроме того, F (.) и H (.) - произвольные линейные или нелинейные функции с компонентами Nx (число состояний) и Ny (количество выходов) соответственно. Можно оценить любой из параметров p1,..., pNpo модели, а также вектор X (0) начального состояния. Стоит подчеркнуть, что
моделирование временных рядов, т.е. моделирование без внешнего входного сигнала u (t), и
статическое моделирование, т.е. моделирование без каких-либо состояний x (t)
это два особых случая, которые поддерживаются IDNLGREY. (Примеры этих двух категорий моделирования см. в учебных пособиях idnlgreydemo3 и idnlgreydemo5).
Первым выполняемым этапом моделирования IDNLGREY всегда является реализация файла модели MATLAB или C MEX, определяющего способ обновления состояний и вычисления выходных данных. Более того, пользователь должен написать файл модели MODFILENAME.m или MODFILENAME.c, определенный со следующими входными и выходными аргументами (обратите внимание, что эта форма требуется для файлов модели типа MATLAB и C MEX).
[dx, y] = MODFILENAME (t, x, u, p1, p2,..., pNpo, FileArgument)
MODFILENAME здесь может быть любым выбранным пользователем именем файла MATLAB или C MEX-файла, например, см. twotanks_m.m, pendulum_c.c и т.д. Этот файл должен быть определен для возврата двух выходных данных:
dx: правая сторона (стороны) уравнения (уравнений) состояния-пространства (вектора столбца с вещественными элементами Nx; [] для статических моделей)
y: правая сторона (стороны) выходного уравнения (уравнений) (вектор столбца с вещественными значениями Ny)
и он должен принимать 3 + Npo (+ 1) входных аргументов, указанных следующим образом:
t: текущее время
x: вектор состояния в момент времени t ([] для статических моделей)
u: входной вектор в момент времени t ([] для моделей временных рядов)
p1, p2,..., pNpo: отдельные параметры (которые могут быть вещественными скалярами, векторами столбцов или 2-мерными матрицами); Npo - это количество объектов параметров, которые для моделей со скалярными параметрами совпадают с числом параметров Np
FileArgument: необязательные входные данные для файла модели
В дальнейшем обсуждении мы сосредоточимся на написании модели с использованием языка MATLAB или файлов C-MEX. Однако IDNLGREY также поддерживает P-файлы (защищенные файлы MATLAB, полученные с помощью команды MATLAB «pcode») и дескрипторы функций. Фактически, возможно использование не только файлов модели C MEX, но и файлов Fortran MEX. Для получения дополнительной информации об этих интерфейсах см. документацию MATLAB по внешним интерфейсам.
Какой файл модели должен быть реализован? Ответ на этот вопрос действительно зависит от использования модели.
Реализация с использованием языка MATLAB (что приводит к созданию файла * .m) имеет определенные преимущества. Во-первых, можно избежать трудоемкого низкоуровневого программирования и сосредоточиться больше на аспектах моделирования. Во-вторых, любая функция, доступная в MATLAB и его панелях инструментов, может использоваться непосредственно в файлах модели. В-третьих, такие файлы будут меньше, и без каких-либо изменений будет автоматически применяться вся встроенная проверка ошибок MATLAB. Кроме того, это получается без какой-либо компиляции кода.
Моделирование C MEX гораздо более активно и требует базовых знаний о языке программирования C. Основным преимуществом файлов модели C MEX является повышенная скорость выполнения. Наш общий совет заключается в проведении моделирования C MEX, когда модель будет использоваться много раз, когда используются большие наборы данных и/или когда структура модели содержит много вычислений. Часто стоит начать с использования файла MATLAB и позже обратиться к аналогу C MEX.
Далее перейдем к моделированию файлов MATLAB и воспользуемся нелинейной структурой модели второго порядка, описывающей систему двух резервуаров, в качестве примера. Подробности моделирования см. в разделе idnlgreydemo2. Содержание twotanks_m.m выглядит следующим образом.
type twotanks_m.mfunction [dx, y] = twotanks_m(t, x, u, A1, k, a1, g, A2, a2, varargin)
%TWOTANKS_M A two tank system.
% Copyright 2005-2006 The MathWorks, Inc.
% Output equation.
y = x(2); % Water level, lower tank.
% State equations.
dx = [1/A1*(k*u(1)-a1*sqrt(2*g*x(1))); ... % Water level, upper tank.
1/A2*(a1*sqrt(2*g*x(1))-a2*sqrt(2*g*x(2))) ... % Water level, lower tank.
];
В заголовке функции мы находим требуемые входные аргументы t, x и u, за которыми следуют шесть параметров скалярной модели, A1, k, a1, g, A2 и a2. В случае файла MATLAB последним входным аргументом всегда должен быть varargin для поддержки передачи необязательного входного аргумента файла модели FileArgument. В объекте модели IDNLGREY FileArgument хранится как массив ячеек, который может содержать любой тип данных. Доступ к первому элементу FileArgument осуществляется через varargin {1} {1}.
Переменные и параметры указываются стандартным способом MATLAB. Первое состояние - x (1) и второе x (2), вход - u (1) (или просто u в случае, если он скалярный), а доступ к скалярным параметрам осуществляется просто через их имена (A1, k, a1, g, A2 и a2). Доступ к отдельным элементам векторных и матричных параметров осуществляют соответственно в виде P (i) (элемента i векторного параметра с именем P) и в виде P (i, j) (элемента в строке i и столбце j матричного параметра с именем P).
Запись файла модели C MEX является более сложной задачей, чем запись файла модели MATLAB. Для упрощения этого шага рекомендуется скопировать доступный шаблон модели IDNLGREY C MEX в MODFILENAME.c. Этот шаблон содержит исходный код скелета, а также подробные инструкции по настройке кода для конкретного приложения. Местоположение файла шаблона определяется путем ввода следующего в командной строке MATLAB.
fullfile(matlabroot, 'toolbox', 'ident', 'nlident', 'IDNLGREY_MODEL_TEMPLATE.c')
В примере с двумя резервуарами этот шаблон был скопирован в twotanks_c.c. После некоторых первоначальных модификаций и конфигураций (описанных ниже) были введены уравнения состояния и вывода, что привело к следующему исходному коду C MEX.
type twotanks_c.c/* Copyright 2005-2015 The MathWorks, Inc. */
/* Written by Peter Lindskog. */
/* Include libraries. */
#include "mex.h"
#include <math.h>
/* Specify the number of outputs here. */
#define NY 1
/* State equations. */
void compute_dx(double *dx, double t, double *x, double *u, double **p,
const mxArray *auxvar)
{
/* Retrieve model parameters. */
double *A1, *k, *a1, *g, *A2, *a2;
A1 = p[0]; /* Upper tank area. */
k = p[1]; /* Pump constant. */
a1 = p[2]; /* Upper tank outlet area. */
g = p[3]; /* Gravity constant. */
A2 = p[4]; /* Lower tank area. */
a2 = p[5]; /* Lower tank outlet area. */
/* x[0]: Water level, upper tank. */
/* x[1]: Water level, lower tank. */
dx[0] = 1/A1[0]*(k[0]*u[0]-a1[0]*sqrt(2*g[0]*x[0]));
dx[1] = 1/A2[0]*(a1[0]*sqrt(2*g[0]*x[0])-a2[0]*sqrt(2*g[0]*x[1]));
}
/* Output equation. */
void compute_y(double *y, double t, double *x, double *u, double **p,
const mxArray *auxvar)
{
/* y[0]: Water level, lower tank. */
y[0] = x[1];
}
/*----------------------------------------------------------------------- *
DO NOT MODIFY THE CODE BELOW UNLESS YOU NEED TO PASS ADDITIONAL
INFORMATION TO COMPUTE_DX AND COMPUTE_Y
To add extra arguments to compute_dx and compute_y (e.g., size
information), modify the definitions above and calls below.
*-----------------------------------------------------------------------*/
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
/* Declaration of input and output arguments. */
double *x, *u, **p, *dx, *y, *t;
int i, np;
size_t nu, nx;
const mxArray *auxvar = NULL; /* Cell array of additional data. */
if (nrhs < 3) {
mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidSyntax",
"At least 3 inputs expected (t, u, x).");
}
/* Determine if auxiliary variables were passed as last input. */
if ((nrhs > 3) && (mxIsCell(prhs[nrhs-1]))) {
/* Auxiliary variables were passed as input. */
auxvar = prhs[nrhs-1];
np = nrhs - 4; /* Number of parameters (could be 0). */
} else {
/* Auxiliary variables were not passed. */
np = nrhs - 3; /* Number of parameters. */
}
/* Determine number of inputs and states. */
nx = mxGetNumberOfElements(prhs[1]); /* Number of states. */
nu = mxGetNumberOfElements(prhs[2]); /* Number of inputs. */
/* Obtain double data pointers from mxArrays. */
t = mxGetPr(prhs[0]); /* Current time value (scalar). */
x = mxGetPr(prhs[1]); /* States at time t. */
u = mxGetPr(prhs[2]); /* Inputs at time t. */
p = mxCalloc(np, sizeof(double*));
for (i = 0; i < np; i++) {
p[i] = mxGetPr(prhs[3+i]); /* Parameter arrays. */
}
/* Create matrix for the return arguments. */
plhs[0] = mxCreateDoubleMatrix(nx, 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(NY, 1, mxREAL);
dx = mxGetPr(plhs[0]); /* State derivative values. */
y = mxGetPr(plhs[1]); /* Output values. */
/*
Call the state and output update functions.
Note: You may also pass other inputs that you might need,
such as number of states (nx) and number of parameters (np).
You may also omit unused inputs (such as auxvar).
For example, you may want to use orders nx and nu, but not time (t)
or auxiliary data (auxvar). You may write these functions as:
compute_dx(dx, nx, nu, x, u, p);
compute_y(y, nx, nu, x, u, p);
*/
/* Call function for state derivative update. */
compute_dx(dx, t[0], x, u, p, auxvar);
/* Call function for output update. */
compute_y(y, t[0], x, u, p, auxvar);
/* Clean up. */
mxFree(p);
}
Рассмотрим содержимое этого файла. В качестве первого наблюдения мы можем разделить работу по написанию файла модели C MEX на четыре отдельных подшаблона, последний из которых является необязательным:
Включение С-библиотек и определение количества выходов.
Запись функции, вычисляющей правую сторону (стороны) уравнения (уравнений) состояния, compute_dx.
Запись функции, вычисляющей правую сторону (стороны) выходного уравнения (уравнений), compute_y.
Опционально обновление функции основного интерфейса, которая включает в себя основные функции проверки ошибок, код для создания и обработки входных и выходных аргументов, а также вызовы compute_dx и compute_y.
Прежде чем мы рассмотрим эти подшаблоны более подробно, кратко прокомментируем пару общих особенностей языка программирования C.
Высокоточные переменные (все входы, состояния, выходы и параметры объекта IDNLGREY) должны иметь тип данных «double».
Унарный оператор *, расположенный непосредственно перед именем переменной или параметра, является так называемым оператором обособления. C-объявление «double * A1;» указывает, что A1 является указателем на двойную переменную. Конструкция указателя - это понятие в C, которое не всегда так легко понять. К счастью, если объявления выходных/входных переменных compute_y и compute_dx не изменены и все распакованные параметры модели внутренне объявлены с помощью *, то нет необходимости знать больше о указателях с точки зрения моделирования IDNLGREY.
compute_y и compute_dx сначала объявляются и реализуются, где после вызываются в функции основного интерфейса. В объявлении ключевое слово «void» явно указывает, что значение не должно быть возвращено.
Дополнительные сведения о языке программирования C см. в книге
B.W. Kernighan and D. Ritchie. The C Programming Language, 2nd
edition, Prentice Hall, 1988.
На первом подэтапе мы сначала включаем C-библиотеки «mex.h» (требуется) и «math.h» (требуется для более продвинутой математики). Количество выходов также объявляется для каждого файла моделирования с использованием стандартного C-define:
/* Include libraries. */
#include "mex.h"
#include "math.h"
/* Specify the number of outputs here. */
#define NY 1
При желании можно также включить большее количество библиотек C, чем библиотеки, описанные выше.
Библиотека «math.h» должна быть включена всякий раз, когда любое состояние или выходное уравнение содержит более продвинутую математику, как тригонометрические и квадратные функции корня. Ниже приведен выбранный список функций, включенных в «math.h» и аналог, найденный в MATLAB:
C-function MATLAB function
========================================
sin, cos, tan sin, cos, tan
asin, acos, atan asin, acos, atan
sinh, cosh, tanh sinh, cosh, tanh
exp, log, log10 exp, log, log10
pow(x, y) x^y
sqrt sqrt
fabs abs
Обратите внимание, что функции MATLAB более универсальны, чем соответствующие C-функции, например, первые обрабатывают комплексные числа, тогда как вторые - нет.
Далее в файле мы найдем функции для обновления состояний, compute_dx и вывода, compute_y. Обе эти функции содержат списки аргументов, при этом выходные данные вычисляются (dx или y) в позиции 1, после чего следуют всем переменным и параметрам, необходимым для вычисления правой стороны (сторон) состояния и выходных уравнений соответственно.
Все параметры содержатся в массиве параметров p. Первый шаг в compute_dx и compute_y состоит в распаковке и присвоении имени параметрам, которые будут использоваться в последующих уравнениях. В twotanks_c.c compute_dx объявляет шесть переменных параметров, значения которых определяются соответствующим образом:
/* Retrieve model parameters. */
double *A1, *k, *a1, *g, *A2, *a2;
A1 = p[0]; /* Upper tank area. */
k = p[1]; /* Pump constant. */
a1 = p[2]; /* Upper tank outlet area. */
g = p[3]; /* Gravity constant. */
A2 = p[4]; /* Lower tank area. */
a2 = p[5]; /* Lower tank outlet area. */
compute_y, с другой стороны, не требует какого-либо параметра для вычисления выходного сигнала, и, следовательно, не извлекается параметр модели.
Как и в случае С, первый элемент массива хранится в позиции 0. Следовательно, dx [0] в C соответствует dx (1) в MATLAB (или просто dx в случае, если это скаляр), вход u [0] соответствует u (или u (1)), параметр A1 [0] соответствует A1 и так далее.
В приведенном выше примере мы используем только скалярные параметры, в этом случае общее количество параметров Np равно количеству объектов параметров Npo. Если какой-либо вектор или параметр матрицы включен в модель, то Npo < Np.
На скалярные параметры ссылаются как P [0] (P (1) или только P в файле MATLAB), а на векторный элемент i: th - как P [i-1] (P (i) в файле MATLAB). Матрицы, передаваемые в файл модели C MEX, отличаются в том смысле, что столбцы расположены друг на друге в очевидном порядке. Следовательно, если P является матрицей 2 на 2, то P (1, 1) упоминается как P [0], P (2, 1) - как P [1], P (1, 2) - как P [2] и P (2, 2) - как P [3]. Пример использования скалярных, векторных и матричных параметров см. в разделе «Учебные пособия по нелинейной идентификации серого ящика: промышленный робот с тремя степенями свободы: C MEX-File Modeling of MIMO System using Vector/Matrix Parameters», idnlgreydemo8.
Функции обновления состояния и вывода могут также включать в себя другие вычисления, чем просто извлечение параметров и вычисление выражений правой стороны. Для скорости выполнения можно, например, объявить и использовать промежуточные переменные, значения которых используются несколько раз в следующих выражениях. Наглядным примером в этом отношении является упомянутое выше учебное пособие по роботам idngreydemo8.
compute_dx и compute_y также могут обрабатывать необязательный FileArgument. Данные FileArgument передаются этим функциям в переменной auxvar, так что первый компонент FileArgument (массив ячеек) может быть получен через
mxArray* auxvar1 = mxGetCell(auxvar, 0);
Здесь mxArray - это определяемый MATLAB тип данных, который обеспечивает обмен данными между C MEX-файлом и MATLAB. В свою очередь, auxvar1 может содержать любые данные. Синтаксический анализ, проверка и использование auxvar1 должны выполняться исключительно в рамках этих функций, где эту функциональность должен реализовать конструктор файлов модели. Для получения дополнительной информации о функциях, работающих с mxArray, рассмотрим документацию MATLAB по внешним интерфейсам. Пример использования дополнительных аргументов файла модели C MEX приведен в разделе idngreydemo6 «Учебные пособия по нелинейной идентификации серого ящика: система передачи сигнала: моделирование файлов C MEX с использованием дополнительных входных аргументов».
Основная функция интерфейса почти всегда должна иметь одинаковое содержимое, и для большинства приложений не требуется никаких изменений. В принципе, единственной частью, которая может быть рассмотрена для изменений, является то, где осуществляются вызовы compute_dx и compute_y. Для статических систем вызов к compute_dx может быть пропущен. В других ситуациях может быть желательным передавать только переменные и параметры, упомянутые в уравнениях состояния и вывода. Например, в уравнении вывода системы двух резервуаров, где используется только одно состояние, можно очень сократить входной список аргументов до
void compute_y(double *y, double *x)
и вызовите compute_y в функции главного интерфейса как
compute_y(y, x);
Списки входных аргументов compute_dx и compute_y также могут быть расширены для включения дополнительных переменных, выведенных в функции интерфейса. Вычисляются и, следовательно, могут передаваться следующие целочисленные переменные: nu (количество входов), nx (количество состояний) и np (здесь количество объектов параметров). Например, nx передается compute_y в модели, исследованной в учебном пособии idnlgreydemo6.
Завершенный файл модели C MEX должен быть скомпилирован, прежде чем его можно будет использовать для моделирования IDNLGREY. Компиляция может быть легко выполнена из командной строки MATLAB как
mex MODFILENAME.c
Обратите внимание, что команда mex-command должна быть настроена до ее первого использования. Это также достигается из командной строки MATLAB через
mex -setup
С готовым к выполнению файлом модели легко создавать объекты модели IDNLGREY, для которых можно выполнять моделирование, оценки параметров и т.д. Мы иллюстрируем это созданием двух различных объектов модели IDNLGREY для описания двух систем резервуаров, один с использованием файла модели, написанного в MATLAB, и один с использованием файла C MEX, подробно описанного выше (обратите внимание, что файл модели C MEX уже скомпилирован).
Order = [1 1 2]; % Model orders [ny nu nx]. Parameters = [0.5; 0.003; 0.019; ... 9.81; 0.25; 0.016]; % Initial parameter vector. InitialStates = [0; 0.1]; % Initial values of initial states. nlgr_m = idnlgrey('twotanks_m', Order, Parameters, InitialStates, 0)
nlgr_m =
Continuous-time nonlinear grey-box model defined by 'twotanks_m' (MATLAB 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 6 free parameter(s) (out of 6).
Status:
Created by direct construction or transformation. Not estimated.
nlgr_cmex = idnlgrey('twotanks_c', Order, Parameters, InitialStates, 0)nlgr_cmex =
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 6 free parameter(s) (out of 6).
Status:
Created by direct construction or transformation. Not estimated.
В этом учебном пособии рассматривается процесс написания файлов моделей IDNLGREY MATLAB и C MEX. Наконец, мы завершаем презентацию перечислением доступных в настоящее время файлов моделей IDNLGREY и учебного пособия/тематического исследования, в которых они используются. Для упрощения дальнейших сравнений мы перечислим как MATLAB (соглашение об именовании FILENAME_m.m), так и файлы модели C MEX (соглашение об именовании FILENAME_c.c) и в столбце учебного пособия укажите тип подхода к моделированию, который используется в учебном пособии или тематическом исследовании.
Tutorial/Case study MATLAB file C MEX-file
======================================================================
idnlgreydemo1 (MATLAB) dcmotor_m.m dcmotor_c.c
idnlgreydemo2 (C MEX) twotanks_m.m twotanks_c.c
idnlgreydemo3 (MATLAB) preys_m.m preys_c.c
(C MEX) predprey1_m.m predprey1_c.c
(C MEX) predprey2_m.m predprey2_c.c
idnlgreydemo4 (MATLAB) narendrali_m.m narendrali_c.c
idnlgreydemo5 (MATLAB) friction_m.m friction_c.c
idnlgreydemo6 (C MEX) signaltransmission_m.m signaltransmission_c.c
idnlgreydemo7 (C MEX) twobodies_m.m twobodies_c.c
idnlgreydemo8 (C MEX) robot_m.m robot_c.c
idnlgreydemo9 (MATLAB) cstr_m.m cstr_c.c
idnlgreydemo10 (MATLAB) pendulum_m.m pendulum_c.c
idnlgreydemo11 (C MEX) vehicle_m.m vehicle_c.c
idnlgreydemo12 (C MEX) aero_m.m aero_c.c
idnlgreydemo13 (C MEX) robotarm_m.m robotarm_c.c
Содержимое этих файлов модели можно отобразить в окне команд MATLAB с помощью команд «type FILENAME_m.m» или «type» FILENAME_c.c. Все файлы модели находятся в папке, возвращенной следующей командой MATLAB.
fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'examples')