Как стать автором
Обновить

Комментарии 13

Не слишком ли велика ошибка численного метода?
Погрешность 28 %. Да, достаточно большая. Причем там появляется систематическая погрешность. Почему именно, нужно разбираться, но мне кажется что эти ошибки возникают из-за использования встроенных алгоритмов MatLab для решения таких обратных задач. По крайней мере оценка хоть примерно сходится в направлении "истинного" значения.
Совершенно не понимаю фразу
примерно на 3,5 секунде процесс стабилизируется
. У вас параметр a оценивается непрерывно во времени? Или когда уже известны обе траектории — модельная и фактическая? Почему такая значительная ошибка в определении а, если шума у вас нет? Самую интересную часть спрятали в код.
Да вы правы, некорректно написал. Параметр оценивается дискретно. И правильнее было сказать, что оценка параметра а сходится к "истинному" значению. Почему такая значительная ошибка? Уже написал выше, я думаю что это из-за использования встроенных алгоритмов для решения обратных задач. Но это мое мнение, если вы считаете по другому, напишите пожалуйста.
У вас нет шума, метод Рунге-Кутты, который вы используете, дает 6-7 верных знаков при стандартных настройках. В моем понимании, оценка для а должна с самого начала совпадать с точным значением (с теми же 6-7 верными знаками). Возможно (пальцем в небо), проблема в том, что вы сравниваете численное решение с аналитическим в разных точках. Метод сравнения мне кажется очень странным, вы сравниваете два решения всего лишь в одной точке (текущей) и по этой точке оцениваете параметр. На мой взгляд, с шумом это будет работать крайне плохо. Более логично сравнивать, например, интеграл квадрата отклонения за весь промежуток времени. Хотя я бы сравнивал локальные куски решения (например, в окрестности текущего момента времени), либо вообще подставлял траекторию в диффур и искал МНК оценку для а.
Точки я использую одни и те же. Эти точки "T" выходят автоматически при решении в самом начале нашей системы методом Рунге-Кутты. Потом я их же использую для параметрической идентификации.
Можете объяснить, почему вы приравниваете производную от x1 по a к значению x1 из модели? У них даже размерность разная. Почему вы игнорируете вторую компоненту решения, хотя она явно фигурирует в МНК?
Потому что сначала мы решаем нашу систему дифф. уравнений

S = dsolve(diff(x) == ax + 1y,'x(0)=20', diff(y) == 0.0225x — 0.3y,'y(0)=20');

И выбираем первую фазовую координату, потому что она содержит наш искомый параметр

x(t) = S.x;

Эта фазовая координата и есть x_1
Затем в цикле

for i=2:81
SSS(i)=solve(SS(i)==X(i,1),a);
end

приравниваем SS (это модельные значения x_1) к X(i,1) (это измеренные значения первой фазовой координаты x_1) для того чтобы найти оценку a^* в моменты времени T

Игнорирую вторую фазовую координату, так как в данный момент она не нужна. И мы можем вообще изначально установить что мы её не наблюдаем.
SS — это производная x по a.
f=diff(x(t),'a');
S1=simplify(f);
SS=eval(S1);

А что мешало, взять оптимальный фильтр Калмана-Бьюси для линейной системы, расширить пространство состояний одним уравнением вашего параметра и решить полученную систему любым методом с любой точностью? и ошибка будет стремится к нулю, для линейной системы точно. пруф — справочник по теории автоматического управления А.А. Красовского.
Я не стремился решить эту задачу наилучшим способом. Просто на этом примере я хотел показать решение задачи параметрической идентификации методом наименьших квадратов , когда модель задана в виде системы линейных дифференциальных уравнений. Безусловно фильтр Калмана лучше справился бы с этой задачей.
Важным свойством оценок по МНК является существование только одного локального минимума, совпадающего с глобальным.

Вообще говоря, это не верно. Минимум будет единственным только в случае выпуклости оптимизируемого функционала. А он будет выпуклым далеко не для всех распределений, поэтому априорные сведения о случайной ошибке все же не будут лишними.

Поэтому оценка a* является единственной.

При фреквентистском подходе к оцениванию, оценка по МНК для выпуклого функционала будет единственной.
Но это, вообще говоря, не единственная оценка искомого параметра a.
А если перейти на байесовы методы, то a вообще является распределением.
Согласен с Вами, спасибо за замечание.
Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.