Численный расчет производной квантованной функции

[_Ivana]

Столкнулся с такой задачей - имеется механическая система, один из параметров которой (угол между двумя звеньями) изменяется относительно плавно, но измеряется датчиком (энкодером), который изменяет значения на выходе ступенчато, с конструктивно заложенным шагом. Количество шагов на полный оборот порядка нескольких тысяч, то есть для измерения угла вполне достаточно. Но требуется также рассчитывать производную этого угла. Если считать как разность текущего и предыдущего показаний деленную на интервал опроса датчика, то получаются весьма большие скачки, которые идут в управление и не способствуют стабильности конструкции. Как можно численно рассчитать производную в точке, имея лишь предыдущие значения квантованной функции?

[TOTAL]

По нескольким последним значениям постройте прямую методом наименьших квадратов.

[_Ivana]

TOTAL спасибо, навскидку это должна быть хорошая идея. Я думал даже про полином 3 степени по МНК по последним точкам, но уже не помню почему не промоделировал этот вариант. Действительно, аппроксимация тут должна гораздо больше подходить, нежели интерполяция, тем более в МНК можно поиграться весами каждой точки. Попробую вечером промоделировать этот вариант.

[Евгений Машеров]

Учитывая то, что распределение ошибки не нормально, а равномерно $U(-\frac h 2, \frac h 2)$, где h - бин датчика, может оказаться лучше не МНК, а чебышевская аппроксимация.

[_Ivana]

Попробовал прямую по МНК, результат вроде получше, чем первая разность на крайних точках из заданного интервала. Насколько я понимаю, и МНК, и первые разности с пропуском любого количество точек, это просто применение к сигналу некоего КИХ фильтра с заданными коэффициентами. Может Чебышев тоже из этой серии. И можно для заданного количества точек подобрать коэффициенты для лучшего результата, поиграться с ними. Но в любом случае, при увеличении временнОго интервала анализа значение производной получается точнее, но у нас появляется задержка - а в моем случае это критично, у меня производная идет прямо в ПИД управление системой, при задержке больше критической система идет вразнос.

Мне рекомендовали еще 2 подхода:
- не равномерный периодический опрос датчика, а максимально точное измерение промежутков времени между щелчками и последующая обработка такого неравномерно дискретизированного сигнала.
- построение математической модели системы и использовать в качестве производной рассчитанные значения из модели, а уже поведение модели корректировать по данным датчика угла.

Оба подхода интересные, но пока я не владею достаточными знаниями для их реализации.

[svv]

_Ivana в сообщении #838471 писал(а):

- не равномерный периодический опрос датчика, а максимально точное измерение промежутков времени между щелчками и последующая обработка такого неравномерно дискретизированного сигнала.

По-моему, самый правильный подход.

[Евгений Машеров]

Если он вообще технически возможен. Обычный АЦП выдаёт равномерный во времени поток отсчётов. Можно, конечно, резко увеличить частоту дискретизации и отслеживать моменты переходов, но тут тоже есть ограничения.
Очень простой вариант - считать производную просто как разности отсчётов, а потом сглаживать экспоненциально
$y_i=\alpha x_i+(1-\alpha)y_{i-1}$

[_Ivana]

У меня контроллер будет работать на частоте нескольких единиц/десятков мегагерц, частота опроса остальных датчиков (гироскопа и акселерометра) - несколько сотен герц. Я могу варьировать частоту опроса остальных датчиков, если за время каждого интервала буду успевать делать нужные вычисления управления в плавучке (или фиксированной точности), и принципиально, что это будет именно периодический опрос. А вот энкодер я повешу на входящее прерывание, и время его щелчков я смогу рассчитывать не с точностью частоты опроса остальных датчиков, а с точностью тактовой частоты МК - а это на много порядков точнее. Я не использую АЦП и не преобразую аналоговый сигнал датчика в цифровой с какой-то частотой дискретизации. У меня энкодер - он сам будет посылать сигналы в точное время своего очередного "щелчка", остается только ловить их. И остается вопрос, как обработать эти данные.
За формулу спасибо, постараюсь понять ее и промоделировать.
И конечно подкупает идея составить математическую модель системы и брать данные из нее, с корректировкой модели по мере поступления щелчков с энкодера.

[_Ivana]

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

ЗЫ если кому интересно, могу показать сравнительные графики.

[svv]

Мне очень интересно, покажите, пожалуйста.
Ещё интересно, какая у Вас угловая скорость (или производная угла). И ограничен ли угол каким-то диапазоном, или эта штука может сколько угодно оборотов сделать.

[_Ivana]

svv, если у вас есть Матлаб, могу выложить код симуляции. Сейчас подготовлю графики, обновлю сообщение. Угол ограничен +- 90 градусами, но не из-за энкодера (тот может обороты крутить бесконечно), просто моя конструкция - балансирующий акробот. Угловая скорость на модели получается довольно большая, боюсь, двигатель такую не обеспечит. Но она зависит от параметров модели. подбором длинн/масс надеюсь вогнать ее в диапазон.

[svv]

_Ivana
А Вы случайно не самокатик типа сегвеевского хотите сделать?
Матлаба, к сожалению, нет. Да я буду рад одной-двум картинкам!

[_Ivana]

Может не совсем информативно, не показано различие в сигнале управления в зависимости от задержки, но сглаживание выбросов и задержка видны явно.
Самокатик сделать имхо проще. Я хочу сделать балансирующий акробот (известный термин), управляемый в межзвенном шарнире. Чтобы поставить его нижним звеном на пол, а он стоял и не падал.

[svv]

Спасибо! А скажите ещё, в каких единицах время на графиках?

[_Ivana]

Я не переходил к безразмерным единицам, параметры модели задавал в СИ (килограммы, метры, секунды). Если я нигде не ошибся, то на графиках время в секундах, а угол в градусах. Графики показывают стабилизацию положения равновесия при начальном отклонении, производная межзвенного угла рассчитывается по 3 вариантам. Есть незначительные отличия в сигнале управления, вызванные как уровнем выбросов, так и задержкой фильтрованной производной.

ЗЫ разумеется, если период опроса угла уменьшить, то на первом графике выбросы будут еще больше, длину КИХ фильтра придется пропорционально увеличить, а БИХ... еще не пробовал, проверю.