В этом примере показано, как предоставить необязательные входные параметры для моделей 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. Следующие три аргумента FileArgument впоследствии будут использоваться при выполнении моделирования 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 переменной auxvar. После объявления и извлечения параметров модели (а также объявления промежуточных переменных, таких как 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.
Под рукой моделируются входно-выходные данные из провода передачи сигнала длиной 1000 м. Эти данные были моделированы с использованием вышеупомянутой агрегированной структуры модели, но использовалась намного большая 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);
Критерий Final Prediction Error (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 на фигура. Благодаря этому модели станут более «стабильными» и симуляции будут «легче».
figure('Name', [z.Name ': step responses']); step(nlgr100, nlgr30, nlgr10); grid on; legend('nlgr100', 'nlgr30', 'nlgr10', 'Location', 'NorthWest');
Фигура 7: Переходные характеристики трех оцененных моделей передачи сигнала IDNLGREY.
Это руководство в первую очередь посвящено использованию необязательных аргументов входного файла при выполнении моделирования IDNLGREY. Мы сосредоточились на типе моделирования C-MEX, но указали, как это сделать и с файлами MATLAB. Для некоторых видов систем оказывается выгодным (повторное использование модели) спроектировать файл модели с большей гибкостью, чем обычно, например, с точки зрения количества состояний, которые нужно использовать. Система передачи сигналов, обработанная в этом руководстве, является одним из примеров в этой категории.