В этом примере показано, как предоставить дополнительные входные параметры моделям IDNLGREY. Обсуждение концентрируется о том, как сделать это для типов C-MEX файлов модели, все же до некоторой незначительной степени, мы также адресуем самые соответствующие параллели к моделированию файла MATLAB.
Основанием для нашего обсуждения будет система передачи сигнала (телеграфный провод) схематично изображенный в следующем рисунке.
Рисунок 1: Схематическое представление системы передачи сигнала.
На расстоянии x (считаемый от того, где входное напряжение применяется) в телеграфном проводе течет во время t ток i (t, x). Соответствующее напряжение является u (t, x) и отношение между текущим, и напряжение может быть описано двумя двойными Дифференциальными уравнениями с частными производными (УЧП):
Уравнения выше часто упоминаются как (вариант) так называемое "Телеграфное уравнение", с L и C быть индуктивностью и длиной емкости на единицу длины, соответственно.
Телеграфное уравнение может быть аппроксимировано системой обыкновенных дифференциальных уравнений, только рассмотрев ток и напряжение в дискретных точках 0, h, 2 расстояния h..., 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, 2 h..., 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) и позвольте выходу быть напряжением в конце провода, i.e., 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 в структуре, и для провода длины 1 000 м, мы попробуем три различных значения 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 реализует модель передачи сигнала. Введите "тип 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 показывают ниже. ("вводят 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 наконец используется для расчета Люфтганза количеств обязательного параметра =-1 / (L*h) и Ch =-1 / (C*h). См. "Примеры на Нелинейной Серой Идентификации Модели Поля: Создание Файлы Модели IDNLGREY" для получения дополнительной информации о FileArgument.
С 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.
Под рукой симулированные данные ввода - вывода из провода передачи сигнала длины 1 000 м. Эти данные были симулированы с помощью вышеупомянутой агрегированной структуры модели, но намного больший 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
Предполагаемые значения параметров не изменяются что как по сравнению с начальными значениями. Как видно ниже, интересный момент - то, что значения параметров двух наименее включенных моделей на самом деле перемещаются еще дальше от истинных значений параметров, тогда как противоположное происходит для наиболее сложной модели. Существует несколько причин того, почему эта сила происходит, e.g., слишком крупная агрегация или что процедура минимизации закончилась со значениями параметров, соответствующими локальному минимуму.
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.
Мы завершаем этот пример путем рассмотрения модульных переходных процессов предполагаемых моделей передачи сигнала. Как видно, все эти модели, кажется, могут взять транспортную задержку. Однако полюса этих линейных моделей все расположены на мнимых осях, i.e., на контуре устойчивости. Это известно, который симуляции таких моделей являются тонкими, чтобы выполнить и могут привести к колеблющимся переходным процессам вида, показанного в этом рисунке.
Основная проблема здесь состоит в том, что этот вариант телеграфного уравнения не моделирует потерь; такие эффекты могут быть включены путем добавления "резисторов" (последовательно с обмотками) и "резисторов проводимости" (параллельно с конденсаторами) к каждому подблоку Люфтганзы-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 также. Для определенных видов систем это оказывается выгодным (повторное использование модели), чтобы спроектировать файл модели с большей гибкостью, чем нормальный, e.g., в терминах количества состояний, чтобы использовать. Система передачи сигнала, обработанная в этом примере, является одним примером в этой категории.