Аппроксимация. Спринт 1. Постановка и изучение задачи.
Работая над проектом контроллера вентиляторов (да, я всё ещё не теряю надежды его закончить, а впоследствии и выпустить 2-ю версию), я наткнулся на яркую и неожиданную задачу, о которой как раз и хочу рассказать в этой статье.
Я подумал что было бы круто отслеживать не только остановку или значительное снижение оборотов управляемых вентиляторов, но и в целом текущее отклонение оборотов от ожидаемых. Т.е. один раз во время сканирования подключенных кулеров снять несколько точек зависимости оборотов от управляющего сигнала, построить по ним аппроксимированную кривую и по ней уже смотреть расхождение текущих показаний с ожидаемыми.
К примеру, есть корпусы ПК, оборудованные передней дверцей, открытие и закрытие которой может оказывать влияние на воздушные потоки, а следовательно и на нагрузку на двигатели. Может быть наоборот, если один вентилятор работает на захват, а другой на выброс воздуха, то они будут облегчать работу друг друга. Теоретически, можно отслеживать затруднения в работе из-за всяких загрязнений (кошачья/собачья шерсть или волосы там намотались, если башня на полу стоит 😬) и всякое такое прочее. Можно и другие бытовые применения придумать. Следить за работой водяной помпы... а можно отслеживать не обороты, а потребляемый ток и вообще любые зависимые друг от друга параметры (при условии, что их зависимость можно описать выбранной аппроксимирующей функцией). Я когда гуглил материалы по теме - встречал на форумах вопросы по аппроксимации показаний датчиков температуры.
OK, c целями разобрались, теперь посмотрим на возможные варианты реализации.
Самый простой и самый бестолковый вариант - аппроксимация прямой по двум точкам. Во-первых, никто нам не гарантирует что именно те две точки, которые мы выберем, будут наиболее подходящими для построения прямой. Поэтому, даже в случае аппроксимации прямой линией гораздо лучше будет выполнить её по множеству точек, а не по двум.
А во-вторых, наша зависимость не обязана быть линейной. Т.е. она может быть линейной, но очень не всегда. Например, если мы управляем двигателем при помощи питания (чистый ШИМ или DC/DC), зависимость оборотов от коэффициента заполнения может быть примерно как на графике ниже:
Поэтому гораздо лучше будет предусмотреть возможность полиномиальной аппроксимации с произвольно задаваемой степенью полинома. Можно выбрать и другие функции на основе которых будет построена аппроксимация, например логарифмическую или экспоненциальную. Но я в качестве базисной функции остановил свой выбор именно на полиноме, который, как мне кажется, для моих задач будет наиболее универсальным вариантом. Ссылки на всякие вспомогательные статьи, калькуляторы и прочий матан я, как обычно, оставлю в конце статьи.
В двух словах, аппроксимация полиномом заключается в подборе такого n и таких коэффициентов a, чтобы функция
принимала значения как можно более близкие к нашим экспериментально полученным точкам. Т.е. в случае с вентиляторами, такая функция должна максимально точно отражать зависимость оборотов от управляющего сигнала, будь это коэффициент заполнения ШИМ или сопротивление дигипота, или выходной сигнал ЦАПа в цепи обратной связи DC/DC-преобразователя.
Очень-очень давно, страшно сказать, больше четверти века назад я изучал численные методы и, как ни странно, до сих пор у меня сохранились смутные воспоминания о предмете. Руководствуясь которыми, я решил для аппроксимации использовать Метод Наименьших Квадратов, а для решения системы линейных уравнений, получающейся на выходе МНК, использовать методы Гаусса и Крамера. Просто потому, что эти методы мне знакомы и они довольно просты для программной реализации... хотя и пришлось немного постараться, чтобы освежить воспоминания.
И у Гаусса и у Крамера есть свои плюсы и свои минусы. В частности, метод Гаусса быстрее и позволяет оперировать меньшими числами, а значит мы с меньшей вероятностью получим значение "inf", если будем использовать float в качестве базового типа. Зато метод Крамера может обойтись вообще без деления, поскольку деление это самая последняя операция в методе. Т.е. можно вообще хранить коэффициенты в виде целочисленных пар числитель/знаменатель и выполнять деление только при вычислении значения полинома в конкретной точке. Красота! Ну и на мой вкус Крамер алгоритмически проще.
Давайте резюмируем задачу: требуется выполнить программную аппроксимацию полиномом n-й степени по m-количеству точек (m и n произвольные) с использованием метода наименьших квадратов и методов Гаусса или Крамера (на выбор). При этом желательно чтобы и расчёт коэффициентов полинома и вычисление аппроксимированных значений могли выполняться средствами не самых мощных микроконтроллеров (вплоть до Mega328P). В общем-то в сети можно найти довольно много готового кода по теме, но хочется, во-первых, как-то в это дело погрузиться и систематизировать. А во-вторых хочется дать работу мозгу, дабы успокоить себя что я не окончательно всё забыл за последние 25+ лет. А бонусом получить собственную библиотеку, учитывающую особенности работы в условиях ограниченных ресурсов и с максимально широкими возможностями по кастомизации.
Что касается численных методов, то сама по себе их реализация особых затруднений не представляет. А метод Гаусса для первичных экспериментов этого спринта я вообще полностью стянул из википедии, благо он на C#, на котором я и планировал для начала потренироваться.
Основных проблем тут две:
1) Выбрать тип данных, в котором будут выполняться вычисления. Обычные целочисленные типы отпадают. Даже если у нас значения x/y исходных точек представлены целыми числами (а это в большинстве случаев будет именно так), то даже расчёт полинома 2-й степени может привести к проблемам с переполнением, не говоря уже о низкой точности целочисленных вычислений. Стало быть, придётся оперировать вещественными типами. Но и тут кроется засада: в Arduino, к примеру, нет типа double, а есть только float. Т.е. сам-то double объявлен, но под ним скрывается всё тот же float. В принципе, для полиномов 2-й или даже 3-й степени этого должно быть достаточно. Будут небольшие проблемки со снижением точности вычислений, но не так чтобы критичные. Но для полиномов 5-й степени float'а скорее всего будет уже маловато. А вот double вполне сгодится. Это значит, что и тип данных для наших расчётов желательно сделать настраиваемым. В простых случаях можно будет обойтись float, получив при этом выигрыш в скорости вычислений и экономию объёма кода и данных. А в тяжёлых случаях выберем double. Тут Arduino автоматически выбывает из соревнования, но для STM32 это не проблема. Немного заглядывая в следующий спринт скажу, что я пробовал на STM32F0 вычислять полином 8-й степени по 52-м точкам на типе double и всё получилось. Сейчас для моих задач мне вряд ли понадобится аппроксимация полиномом со степенью больше тройки. Ну, максимум 5-я степень. Но на будущее всегда полезно иметь запас.
И если уж говорить про "запас", то есть ещё один интересный вариант: использовать рациональный тип данных, т.е. дробь с целыми числителем и знаменателем. И все вычисления выполнять именно с дробями. Естественно, что в этом случае тоже в полный рост встаёт проблема с переполнением: обычные целочисленные типы таких вычислений сделать не позволят даже если сокращать дроби после (или даже во время) каждого действия. Решить эту проблему можно используя "длинную арифметику", когда числа представляются не одним значением фиксированной ширины (byte, short, int, long), а произвольным набором знаков, так называемых "лимбов". Т.е. каждое число фактически будет представлено вектором. Вычисления в этом случае выполняются самым обычным школьным "столбиком". Но не все. Поделить в столбик вектор можно только на 1 лимб. А в случае с делением вектора на вектор начинается отдельный ад, трэш и содом. Но об этом более подробно будет в следующем спринте)))
В общем, действительно очень любопытный вариант. Основной и бесспорный плюс которого - идеальность вычислений. Никаких потерь и неточностей (за исключением того, что сам по себе метод наименьших квадратов не идеален). Это ли не круто?!!
Для длинной арифметики существует несколько библиотек, самая известная из которых это GMP (The GNU MP Bignum Library), а точнее, входящая в GMP компактная библиотека mini-GMP. Но для наших экспериментов есть вариант гораздо лучше - воспользоваться готовой C#-библиотекой Fractions, основанной на типе BigInteger - реализации длинной арифметики от Microsoft (исходный код открыт). Это сильно упрощает погружение в тему и, конечно, именно этим путём мы и пойдём! Я понимаю что скорее всего применение дробей и длинных вычислений излишне, но уж очень хочется заиметь себе супер-компактный вариант mini-GMP для микроконтроллеров, этакий tiny-GMP.
2) Я предполагаю, что ограниченность ресурсов МК главным образом выразится в нехватке оперативной памяти. По большому счёту, на время работы алгоритма вычисления коэффициентов мне наплевать. Эта операция должна производиться только один раз, причём либо во время первичной инициализации какого-либо устройства (подключили моторчик к ардуине, получили данные по 10-20 точкам, рассчитали коэффициенты и сохранили в EEPROM), либо в худшем случае эта операция будет выполняться сразу после включения устройства и потратить на неё 2-3-5 секунд - вообще не проблема, особенно по сравнению со временем, которое придётся потратить на получение первичных показаний. Хотя, оптимизировать скорость работы всё же желательно. А вот объём занимаемой памяти это как раз проблема. Алгоритм во время работы должен будет хранить матрицу размером n+2 на n+1 и это надо учитывать. А особенно это придётся учитывать для длинной арифметики. Дело в том, что подобные проекты основаны на технике динамического выделения памяти под число. Допустим, выполнили мы операцию 982359872341 *= 34917540 и в результате неё первое число увеличилось на какое-то количество знаков. При этом новый вектор под большее количество знаков выделяется при помощи realloc(). Ну и вообще, под всеми операциями с длинными числами лежат (в самом простом случае) malloc()/realloc()/free(). А это значит, что даже если у нас мощный МК, на котором мы можем не особенно напрягаясь скомпилировать оригинальную mini-GMP, то в результате большого количества операций мы неизбежно получим сильную фрагментацию памяти. На ESP32 и старших сериях STM32 это вполне себе прокатит, но для более скромных MCU лучше так не делать.
Поэтому, в качестве рабочего решения, я пока думаю выделять фиксированный объём памяти под длинные числа и контролировать переполнение. Посмотрим что получится... Пока могу только сказать что предварительные эксперименты по следующему спринту показали, что базового типа, представленного массивом из 4-х uint32_t (16 байт на число) вполне хватает для вычисления коэф-тов полинома 3-й степени. Но это сильно зависит и от исходных данных и от метода (с Крамером хуже, но и Гаусс не гарантирует, что дроби будут сокращаться). В общем, тема интересная и я постараюсь за следующий спринт накидать какую-то библиотеку на C++ и оценить насколько решение с дробями и длинной арифметикой будет применимо для эмбедированных систем среднего и младшего уровней.
Но в любом случае, где-нибудь и когда-нибудь мне это пригодится, а может быть и кому-то из вас тоже, так что я считаю что время будет потрачено не зря.
Как итог этого спринта я прикладываю исходный код консольной C#-программы, которая по заданным экспериментальным точкам (там 3 варианта, но вы можете попробовать любые свои) вычисляет коэффициенты полинома заданной степени сразу обоими методами и приводит оценку корреляции и процента ошибки для каждого из них.
Базовый тип данных для вычислений можно выбрать любой из трёх: Single (аналог float в C#), Double или Fraction (основанный на BigInteger). Получившийся полином вполне можно подставить в Desmos или Excel и визуально оценить совпадение экспериментальных точек с аппроксимирующей кривой.
Поскольку это C#, то он легко компилируется и на Линукс и на Маке и, само-собой на ПеКа, так что с этим проблем быть не должно. Только не забудьте поставить библиотеку Fractions, ссылка на неё в самом низу.
На этом первый вводный спринт я могу закрыть. Получилось не очень объёмно, но на то он и вводный спринт. Буду держать вас в курсе.
Метод наименьших квадратов - http://vmath.ru/vf5/interpolation/mnk
Метод наименьших квадратов - https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%BD%D0%B0%D0%B8%D0%BC%D0%B5%D0%BD%D1%8C%D1%88%D0%B8%D1%85_%D0%BA%D0%B2%D0%B0%D0%B4%D1%80%D0%B0%D1%82%D0%BE%D0%B2
Метод Крамера - https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%9A%D1%80%D0%B0%D0%BC%D0%B5%D1%80%D0%B0
Метод Гаусса - https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%93%D0%B0%D1%83%D1%81%D1%81%D0%B0
Определитель матрицы - https://ru.wikipedia.org/wiki/%D0%9E%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B8%D1%82%D0%B5%D0%BB%D1%8C
Калькулятор МНК - https://planetcalc.ru/8735/
Калькулятор Метод Крамера - https://planetcalc.ru/6007/
Калькулятор Метод Гаусса - https://planetcalc.ru/3571/
Калькулятор целочисленный Метод Гаусса (вычисления в форме дробей) - https://planetcalc.ru/6669/
Пример класса для работы с дробями - https://teletype.in/@selfinstallation/3GBKzLBEB
Тип данных C# BigInteger - https://github.com/microsoft/referencesource/blob/master/System.Numerics/System/Numerics/BigInteger.cs
Библиотека C# Fractions - https://www.nuget.org/packages/fractions/