В этом примере показано, как предоставить необязательные входные аргументы моделям IDNLGREY. Обсуждение сосредоточено на том, как это сделать для файлов моделей C-MEX, но в некоторой степени мы также рассмотрим наиболее актуальные параллели с моделированием файлов MATLAB.
Основой для нашего обсуждения будет система передачи сигнала (телеграфный провод), схематически изображенная на следующем рисунке.
Рис. 1. Схематический вид системы передачи сигналов.
На расстоянии x (отсчитываемом от места приложения входного напряжения) в телеграфном проводе протекает в момент t тока i (t, x). Соответствующее напряжение равно u (t, x), и зависимость между током и напряжением может быть описана двумя связанными дифференциальными уравнениями в частных производных (PDE):
Вышеприведенные уравнения часто упоминаются как (вариант) так называемого «телеграфного уравнения», при этом L и C являются индуктивностью и емкостью на единицу длины соответственно.
Телеграфное уравнение может быть аппроксимировано системой обыкновенных дифференциальных уравнений только с учетом тока и напряжения в дискретных точках расстояния 0, h, 2h,..., Nh, где h - расстояние дискретизации и N - общее число таких расстояний. После этого приближения можно считать, что провод состоит из ряда структурно равных сегментов, соединенных друг с другом в цепочку. В литературе этот тип аппроксимации обычно называют агрегацией.
Пусть напряжение и ток на расстоянии x = kh, для k = 0, 1,..., N, в момент времени t обозначаются u_k (t) и i_k (t) соответственно. Аппроксимировать d/dx u (t, x) и d/dx i (t, x) через следующие простые аппроксимации разности:
d u(t, x) u_{k+1}(t) - u_k(t)
--------- ~ ------------------- Forward approximation.
dx h d i(t, x) i_k(t) - i_{k-1}(t)
--------- ~ ------------------- Backward approximation.
dx hдля x = kh.
Причина использования прямого приближения для d/dx u (t, x) заключается в том, что i_N (t) = 0 (разомкнутый провод), так что оставшиеся N дискретизированных токов могут быть смоделированы следующими N дифференциальными уравнениями:
d i_k(t) 1
-------- = - -- (u_{k+1}(t) - u_k(t)) for k = 0, 1, 2, ..., N-1
dt LhАналогично, поскольку u_0 (t) является известным входным сигналом и для его описания не требуется дифференциального уравнения, удобно использовать схему обратного приближения для моделирования d/dx i (t, x) в точках h, 2h,..., N:
d u_k(t) 1
-------- = - -- (i_k(t) - i_{k-1}(t)) for k = 1, 2, ..., N-1
dt Ch d u_N(t) 1 1
-------- = - -- (i_N(t) - i_{N-1}(t) = -- i_{N-1}(t)
dt Ch ChТаким образом, мы пришли к модельной аппроксимации, показанной на фиг.1, в которой уравнения выражены в виде ряда взаимосвязанных катушек и конденсаторов.
Представим 2 * N состояний x1 (t) = i_0 (t), x2 (t) = u_1 (t), x3 (t) = i_1 (t), x4 (t) = u_2 (t),..., x2N-1 (t) = i_{N-1} (t), x2N (t) = u_N (t). Также обозначим вход u (t) = u_0 (t) и пусть выходом будет напряжение на конце провода, т.е. y (t) = x2N (t) = u_N (t). Кажущиеся замены, наконец, приводят к следующей структуре модели пространства состояния.
x1(t) = -1/(Lh)*(x2(t) - u(t))
x2(t) = -1/(Ch)*(x3(t) - x1(t))
x3(t) = -1/(Lh)*(x4(t) - x2(t))
x4(t) = -1/(Ch)*(x5(t) - x3(t))
...
x2N-1(t) = -1/(Lh)*(x2N(t) - x2N-2(t))
x2N(t) = 1/(Ch)*x2N-1(t)y(t) = x2N(t)
В целом, вышеуказанные математические манипуляции привели нас к стандартной структуре модели пространства линейного состояния, с которой вполне может справиться IDGREY - линейный аналог IDNLGREY. Мы не будем выполнять никакое моделирование IDGREY здесь, но вместо этого покажем, как дополнительные входные аргументы могут использоваться IDNLGREY для повышения гибкости моделирования.
Для получения гибкости здесь желательно иметь объект модели IDNLGREY, который сразу же способен иметь дело с проводом любой длины L. Для целей качества моделирования также должно быть простым изменять количество агрегированных блоков N так, чтобы можно было получить достаточно хорошую системную аппроксимацию. Эти требования могут быть обработаны посредством передачи N и L в свойстве FileArgument объекта IDNLGREY. Свойство FileArgument должно быть массивом ячеек, но этот массив может содержать любой тип данных. В этой заявке мы выбираем, чтобы обеспечить N и L в структуре, и для провода длиной 1000 м мы попробуем три различных значения N: 10, 30 и 100. Следующие три параметра FileArguments будут впоследствии использоваться при выполнении моделирования IDNLGREY:
FileArgument10 = {struct('N', 10, 'L', 1000)}; % N = 10 --> 20 states.
FileArgument30 = {struct('N', 30, 'L', 1000)}; % N = 30 --> 60 states.
FileArgument100 = {struct('N', 100, 'L', 1000)}; % N = 100 --> 200 states.
Синтаксический анализ и проверка данных, содержащихся в FileArgument, должны выполняться в файле модели IDNLGREY, где эту функциональность должен реализовать конструктор файлов модели. Для получения h и N в функции без проверки ошибок можно использовать следующие команды:
N = varargin{1}{1}.N;
h = varargin{1}{1}.L/N;Обратите внимание, что параметр FileArgument здесь соответствует параметру varargin {1}, который является последним аргументом, переданным в файл модели. Файловый signaltransmission_m.m реализует модель передачи сигнала. Введите «type signaltransmission_m.m» для просмотра всего файла.
Ситуация несколько более интересна при использовании файлов модели C-MEX, как будет сделано дальше. В этом случае мы заранее не знаем количество состояний Nx, но оно вычисляется в функции главного интерфейса и, таким образом, может быть передано в качестве входного аргумента compute_dx и compute_y. Описание этих функций становится следующим:
void compute_dx(double *dx, int nx, double *x, double *u, double **p, const mxArray *auxvar) void compute_y(double *y, int nx, double *x)
где мы забрали t-переменную (не используется в уравнениях) и вместо этого включили целое число nx в качестве второго аргумента для обеих этих функций. Кроме того, мы также удалили стандартные три конечных аргумента compute_y, поскольку они не используются для вычисления выходных данных. При этих изменениях compute_dx и compute_y вызываются из функции главного интерфейса следующим образом.
/* Call function for state derivative update. */ compute_dx(dx, nx, x, u, p, auxvar);
/* Call function for output update. */ compute_y(y, nx, x);
Первая часть compute_dx показана ниже. («type signaltransmission_c.c» отображает весь файл в окне команд MATLAB ®.)
/* State equations. */
void compute_dx(double *dx, int nx, double *x, double *u, double **p,
const mxArray *auxvar)
{
/* Declaration of model parameters and intermediate variables. */
double *L, *C; /* Model parameters. */
double h, Lh, Ch; /* Intermediate variables/parameters. */
int j; /* Equation counter. */ /* Retrieve model parameters. */
L = p[0]; /* Inductance per unit length. */
C = p[1]; /* Capacitance per unit length. */ /* Get and check FileArgument (auxvar). */
if (mxGetNumberOfElements(auxvar) < 1) {
mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidFileArgument",
"FileArgument should at least hold one element.");
} else if (mxIsStruct(mxGetCell(auxvar, 0)) == false) {
mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidFileArgument",
"FileArgument should contain a structure.");
} else if ( (mxGetFieldNumber(mxGetCell(auxvar, 0), "N") < 0)
|| (mxGetFieldNumber(mxGetCell(auxvar, 0), "L") < 0)) {
mexErrMsgIdAndTxt("IDNLGREY:ODE_FILE:InvalidFileArgument",
"FileArgument should contain a structure with fields 'N' and 'L'.");
} else {
/* Skip further error checking to obtain execution speed. */
h = *mxGetPr(mxGetFieldByNumber(mxGetCell(auxvar, 0), 0,
mxGetFieldNumber(mxGetCell(auxvar, 0), "L"))) / (0.5*((double) nx));
}
Lh = -1.0/(L[0]*h);
Ch = -1.0/(C[0]*h);
...Стоит подчеркнуть, что FileArgument передается compute_dx в переменном ауксваре. После объявления и извлечения параметров модели (а также объявления промежуточных переменных, таких как h, auxvar проверяется на непротиворечивость через ряд внешних интерфейсных подпрограмм (так называемых mx-подпрограмм). Эти подпрограммы MATLAB позволяют создавать, получать доступ, управлять и уничтожать переменные mxArray. Для получения дополнительной информации см. документацию MATLAB по внешним интерфейсам. Последнее предложение else выполняется, если все проверки пройдены, и в этом случае значение поля auxvar N используется для определения значения h. Он и параметры модели L и C в конечном итоге используются для вычисления требуемых величин параметров Lh = -1/( L * h) и Ch = -1/( C * h). Дополнительные сведения о FileArgument см. в разделе «Учебные пособия по идентификации модели нелинейного серого ящика: создание файлов модели IDNLGREY».
При вычисленных Lh и Ch вторая часть compute_dx выглядит следующим образом. Обратите внимание, в частности, на то, как nx используется в цикле for compute_dx для определения количества состояний текущей модели.
...
/* x[0] : Current i_0(t). */
/* x[1] : Voltage u_1(t). */
/* x[2] : Current i_1(t). */
/* x[3] : Voltage u_1(t). */
/* ... */
/* x[Nx-2]: Current i_Nx-1(t). */
/* x[Nx-1]: Voltage u_Nx(t). */
for (j = 0; j < nx; j = j+2) {
if (j == 0) {
/* First transmitter section. */
dx[j] = Lh*(x[j+1]-u[0]);
dx[j+1] = Ch*(x[j+2]-x[j]);
} else if (j < nx-3) {
/* Intermediate transmitter sections. */
dx[j] = Lh*(x[j+1]-x[j-1]);
dx[j+1] = Ch*(x[j+2]-x[j]);
} else {
/* Last transmitter section. */
dx[j] = Lh*(x[j+1]-x[j-1]);
dx[j+1] = -Ch*x[j];
}
}
}Выходная функция обновления compute_dy намного проще функции обновления состояния:
/* Output equation. */
void compute_y(double *y, int nx, double *x)
{
/* y[0]: Voltage at the end of the transmitter. */
y[0] = x[nx-1];
}Теперь мы находимся в положении, в котором мы прямо можем ввести вышеуказанную информацию в три различных объекта IDNLGREY, один для N = 10, один для N = 30 и один для N = 100. Обратите внимание, что отличиями при создании этих моделей являются порядок, вектор начального состояния и используемый необязательный входной аргумент, но именно он.
FileName = 'signaltransmission_c'; % File describing the model structure. Parameters = struct('Name', {'Inductance per unit length' ... % Initial parameters. 'Capacitance per unit length'}, ... 'Unit', {'H/m' 'F/m'}, ... 'Value', {0.99e-3 0.99e-3}, ... 'Minimum', {eps(0) eps(0)}, ... % L, C > 0! 'Maximum', {Inf Inf}, ... 'Fixed', {false false}); % A. Signal transmission model with N = 10; Order10 = [1 1 2*FileArgument10{1}.N]; % Model orders [ny nu nx]. InitialStates10 = zeros(2*FileArgument10{1}.N, 1); % Initial intitial states. nlgr10 = idnlgrey(FileName, Order10, Parameters, InitialStates10, 0, ... 'FileArgument', FileArgument10, ... 'Name', '10 blocks', 'TimeUnit', 's'); % B. Signal transmission model with N = 30; Order30 = [1 1 2*FileArgument30{1}.N]; % Model orders [ny nu nx]. InitialStates30 = zeros(2*FileArgument30{1}.N, 1); % Initial value of intitial states. nlgr30 = idnlgrey(FileName, Order30, Parameters, InitialStates30, 0, ... 'FileArgument', FileArgument30, ... 'Name', '30 blocks', 'TimeUnit', 's'); % C. Signal transmission model with N = 100; Order100 = [1 1 2*FileArgument100{1}.N]; % Model orders [ny nu nx]. InitialStates100 = zeros(2*FileArgument100{1}.N, 1); % Initial value of intitial states. nlgr100 = idnlgrey(FileName, Order100, Parameters, InitialStates100, 0, ... 'FileArgument', FileArgument100, ... 'Name', '100 blocks', 'TimeUnit', 's');
Имена и единицы трех моделей передачи сигнала также устанавливаются, после чего размеры моделей текстуально подтверждаются. Вход - это напряжение, приложенное к проводу, а выход - это напряжение на конце провода.
set(nlgr10, 'InputName', 'Vin', 'InputUnit', 'V', ... 'OutputName', 'Vout', 'OutputUnit', 'V'); set(nlgr30, 'InputName', nlgr10.InputName, 'InputUnit', nlgr10.InputUnit, ... 'OutputName', nlgr10.OutputName, 'OutputUnit', nlgr10.OutputUnit); set(nlgr100, 'InputName', nlgr10.InputName, 'InputUnit', nlgr10.InputUnit, ... 'OutputName', nlgr10.OutputName, 'OutputUnit', nlgr10.OutputUnit); nlgr10 nlgr30 nlgr100
nlgr10 =
Continuous-time nonlinear grey-box model defined by 'signaltransmission_c' (MEX-file):
dx/dt = F(t, u(t), x(t), p1, p2, FileArgument)
y(t) = H(t, u(t), x(t), p1, p2, FileArgument) + e(t)
with 1 input(s), 20 state(s), 1 output(s), and 2 free parameter(s) (out of 2).
Name: 10 blocks
Status:
Created by direct construction or transformation. Not estimated.
nlgr30 =
Continuous-time nonlinear grey-box model defined by 'signaltransmission_c' (MEX-file):
dx/dt = F(t, u(t), x(t), p1, p2, FileArgument)
y(t) = H(t, u(t), x(t), p1, p2, FileArgument) + e(t)
with 1 input(s), 60 state(s), 1 output(s), and 2 free parameter(s) (out of 2).
Name: 30 blocks
Status:
Created by direct construction or transformation. Not estimated.
nlgr100 =
Continuous-time nonlinear grey-box model defined by 'signaltransmission_c' (MEX-file):
dx/dt = F(t, u(t), x(t), p1, p2, FileArgument)
y(t) = H(t, u(t), x(t), p1, p2, FileArgument) + e(t)
with 1 input(s), 200 state(s), 1 output(s), and 2 free parameter(s) (out of 2).
Name: 100 blocks
Status:
Created by direct construction or transformation. Not estimated.
Эти данные были смоделированы с использованием вышеуказанной агрегированной структуры модели, но было использовано значительно большее число N (1500). Моделирование проводили в течение 20 секунд, используя частоту отбора проб 0,1 секунды. Используемые параметры модели были L = C = 1e-3, а начальной точкой был провод нулевого напряжения (все начальные состояния, таким образом, равны нулю). Оба параметра модели, по общему признанию, намного выше, чем для типичного провода передачи сигнала, но были выбраны таким образом, чтобы лучше проиллюстрировать задержку передачи для систем этого типа. Мы загружаем эти данные и помещаем их в объект IDDATA z:
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'signaltransmissiondata')); z = iddata(vout, vin, 0.1, 'Name', 'Signal transmission', ... 'InputName', 'Vin', 'InputUnit', 'V', ... 'OutputName', 'Vout', ... 'OutputUnit', 'V', 'Tstart', 0, 'TimeUnit', 's');
График входных-выходных данных ясно показывает задержку переноса от приложенного напряжения к выходному напряжению на конце провода. Первый импульс входного напряжения на расстоянии около 1,4 секунды появляется на выходе примерно через секунду.
figure('Name', [z.Name ': input-output data']); plot(z);
Рис. 2: Данные ввода-вывода из системы передачи сигнала.
Насколько хорошо работают эти три первоначальные модели? Давайте изучим это, проводя моделирование модели с использованием COMPARE. С точки зрения исполнения, здесь жизненно важно, чтобы начальные состояния не оценивались, таким образом, выбор нулевых векторов начального состояния. Причина этого в том, что число состояний очень велико, особенно для nlgr100 (= 200), и что оценка векторов начального состояния приведет к действительно длительным вычислениям.
compare(z, nlgr100, nlgr30, nlgr10, compareOptions('InitialCondition', 'zero'));
Рис. 3: Сравнение истинного выхода с моделируемыми выходами трех исходных моделей передачи сигнала.
С точки зрения соответствия существует значительная разница между тремя моделями, и, как и ожидалось, наиболее сложная модель превосходит две другие. Разницу в мощности моделирования, возможно, лучше всего рассматривать, глядя на ошибки прогнозирования каждой модели.
peOpt = peOptions('InitialCondition','zero'); e = {pe(z, nlgr100, peOpt) pe(z, nlgr30, peOpt) pe(z, nlgr10, peOpt)}; figtitle = {'nlgr100 : e@' 'nlgr30 : e@' 'nlgr10 : e@'}; figure('Name', [z.Name ': prediction errors']); for i = 1:3 subplot(3, 1, i); plot(e{i}.SamplingInstants, e{i}.OutputData, 'r'); title(['Initial ' figtitle{i} z.OutputName{1}]); axis('tight'); end
Рис. 4: Ошибки прогнозирования, полученные в трех исходных моделях передачи сигналов IDNLGREY.
Хотя здесь мы начали с почти правильных значений параметров истинной системы, это не обязательно означает, что эти значения являются теми, которые обеспечивают наилучшее соответствие модели. Рассмотрим это, оценив два параметра модели трех моделей передачи сигнала с помощью NLGREYEST. Эти вычисления могут занять некоторое время.
opt = nlgreyestOptions('Display', 'on', 'EstimateCovariance', false); nlgr100 = nlgreyest(z, nlgr100, opt); nlgr30 = nlgreyest(z, nlgr30, opt); nlgr10 = nlgreyest(z, nlgr10, opt);
Критерий ошибки окончательного прогнозирования (FPE), применяемый к оценочным моделям, продолжает указывать на то, что наиболее сложная модель превосходит:
fpe(nlgr100, nlgr30, nlgr10)
0.5228 0.6771 4.7749
Оценочные значения параметров изменяются не так сильно по сравнению с начальными значениями. Как можно видеть ниже, интересным моментом является то, что значения параметров двух наименее задействованных моделей фактически отходят дальше от истинных значений параметров, тогда как для наиболее сложной модели происходит обратное. Существует несколько причин, по которым это может произойти, например, слишком грубая агрегация или то, что процедура минимизации завершилась значениями параметров, соответствующими локальному минимуму.
disp(' True nlgr100 nlgr30 nlgr10'); ptrue = [1e-3; 1e-3]; fprintf(' %1.7f %1.7f %1.7f %1.7f\n', ... [ptrue'; getpvec(nlgr100)'; getpvec(nlgr30)'; getpvec(nlgr10)']);
True nlgr100 nlgr30 nlgr10 0.0010000 0.0009952 0.0009759 0.0009879 0.0010000 0.0009952 0.0009759 0.0009879
Затем мы снова используем COMPARE для выполнения моделирования трех оцененных моделей передачи сигнала.
figure compare(z, nlgr100, nlgr30, nlgr10, compareOptions('InitialCondition', 'zero'));
Рис. 5: Сравнение истинного выходного сигнала с моделируемыми выходными сигналами трех расчетных моделей передачи сигнала.
Во всех трех случаях мы получаем немного улучшенную посадку. Хотя это улучшение невелико, оно может быть распознано путем сравнения ошибок прогнозирования, полученных после и до оценки параметров.
peOpt = peOptions('InitialCondition','zero'); e = {pe(z, nlgr100, peOpt) pe(z, nlgr30, peOpt) pe(z, nlgr10, peOpt)}; figtitle = {'nlgr100 : e@' 'nlgr30 : e@' 'nlgr10 : e@'}; figure('Name', [z.Name ': prediction errors']); for i = 1:3 subplot(3, 1, i); plot(e{i}.SamplingInstants, e{i}.OutputData, 'r'); title(['Estimated ' figtitle{i} z.OutputName{1}]); axis('tight'); end
Рис. 6: Ошибки прогнозирования, полученные с помощью трех расчетных моделей передачи сигналов IDNLGREY.
Мы завершаем этот учебник изучением единичных ступенчатых откликов расчетных моделей передачи сигнала. Как видно, все эти модели, по-видимому, способны уловить задержку транспортировки. Однако все полюса этих линейных моделей расположены на воображаемых осях, т.е. на границе устойчивости. Хорошо известно, что моделирование таких моделей является деликатным для выполнения и может привести к колебательным ступенчатым откликам типа, показанного на этой фигуре.
Главная проблема здесь в том, что этот вариант телеграфного уравнения не моделирует никаких потерь; такие эффекты могут быть включены добавлением «резисторов» (последовательно с катушками) и «резисторов проводимости» (параллельно с конденсаторами) к каждому Lh-Ch подблоку фигуры 1. К этому модели станут более «стабильными» и моделирование будет «проще».
figure('Name', [z.Name ': step responses']); step(nlgr100, nlgr30, nlgr10); grid on; legend('nlgr100', 'nlgr30', 'nlgr10', 'Location', 'NorthWest');
Рис. 7: Пошаговые ответы трех оцененных моделей передачи сигналов IDNLGREY.
В данном учебном пособии в первую очередь рассматривается использование дополнительных аргументов входных файлов при моделировании IDNLGREY. Мы сосредоточились на моделировании типа C-MEX, но указали, как это сделать с файлами MATLAB. Для определенных типов систем оказывается выгодным (повторное использование модели) конструировать файл модели с большей гибкостью, чем обычно, например, с точки зрения количества используемых состояний. Система передачи сигналов, рассматриваемая в данном учебном пособии, является одним из примеров в этой категории.