Хотите преодолеть разрыв между математической моделью, которая у вас есть на бумаге, и вашими реальными данными об оборудовании? Python упрощает калибровку модели!

Представьте себе: я даю вам случайную механическую пружину и прошу дать мне математическую модель, описывающую поведение пружины. То есть я хочу, чтобы вы дали мне уравнение, которое скажет мне, какое усилие будет производить пружина при определенном растяжении пружины. Из ваших предыдущих исследований вы знаете, что пружины обычно моделируются в соответствии с законом Гука следующим образом:

где F — выходная сила пружины, k — жесткость пружины, а x — величина растяжения пружины.

Чтобы полностью описать поведение пружины, вам просто нужно определить жесткость пружины k. Итак, вы разрабатываете эксперимент (используя устройство, подобное показанному ниже), чтобы измерить и записать силу, выдаваемую пружиной при различных шагах растяжения x.

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

Так почему я заставил вас пробежаться по этому простому примеру? Потому что он иллюстрирует более общую проблему, которая часто возникает в инженерии — калибровку модели. То есть, учитывая некоторый набор аппаратных средств (в данном случае пружину) и предлагаемую математическую модель этого аппаратного обеспечения (F = kx) с одним или несколькими настраиваемыми параметрами (k), найдите наилучший набор настраиваемых параметров, чтобы модель точно воспроизводит реакцию оборудования. По сути, вы калибруете/настраиваете математическую модель так, чтобы модель и оборудование вели себя одинаково.

Но почему нас это волнует? Математические модели служат основой, на которой мы выполняем нашу инженерную работу! Они дают нам представление о том, какие параметры мы можем спроектировать или изменить, позволяют моделировать систему без проведения дорогостоящих или, возможно, опасных экспериментов, а также позволяют систематически разрабатывать сложные системы управления. Несмотря на то, что приведенный выше пример относительно прост, калибровка модели, как правило, непростая задача, и процедура, необходимая для определения параметров модели, обычно варьируется от системы к системе. Несмотря на это, язык программирования Python почти наверняка можно использовать для упрощения процесса.

В этой статье мы рассмотрим подробный пример того, как можно использовать Python для выполнения калибровки модели в немного более сложной системе. Представленный пример затронет широкий спектр навыков Python, включая объектно-ориентированное программирование, интерполяцию данных, системное моделирование и оптимизацию. Хотя приведенный пример охватывает калибровку модели резистивно-индуктивной (RL) нагрузки, обсуждаемые методологии и методы кодирования широко применимы во всех областях техники.

Проблема

Допустим, у нас есть RL-нагрузка, подобная той, что показана ниже — она может представлять собой электродвигатель или какой-либо другой привод, поведение которого мы хотели бы смоделировать. Нам нужно определить номинал резистора (R) и катушки индуктивности (L). Самый простой способ сделать это — просто подключить LCR-метр и измерить параметры схемы. Но что, если у вас его нет? И что вообще делает счетчик?

Стандартный метод калибровки модели

LCR-метр не творит никакой магии — под капотом он просто подает возбуждение по напряжению v(t), измеряет текущую характеристику i(t), а затем возвращает параметры, которые должны быть для получения этой характеристики. Это та же идея, что и в приведенном выше примере пружины — применить возбуждение (силу), измерить отклик (смещение), а затем вернуть параметры модели (жесткость). Что делает эту проблему более сложной, чем проблема с пружиной, так это то, что теперь реакция системы зависит от времени (подача постоянного напряжения будет производить ток, изменяющийся во времени). Причина этого становится понятной, когда мы формулируем уравнения, описывающие поведение схемы. Применяя законы Кирхгофа, получаем следующее уравнение:

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

Однако можно сделать ключевое наблюдение: как только форма сигнала входного напряжения определена, решение для тока становится алгебраическим выражением. Основываясь на этом наблюдении, один из методов калибровки нашей модели может быть следующим: 1) Выбрать некоторую общую функцию возбуждения для напряжения (возможно, ступенчатую функцию или синусоиду) с известным решением для тока. 2) Подайте это напряжение на цепь и измерьте отклик по току. 3) Сопоставьте кривую аналитического уравнения токового отклика с данными, чтобы исключить параметры схемы. Эта процедура продемонстрирована на графике ниже для фиктивной нагрузки, к которой приложено ступенчатое напряжение — подгоняя экспоненциальную кривую к текущим данным, мы можем оценить индуктивность и сопротивление.

Более общий метод калибровки модели

Вышеупомянутая процедура полностью действительна и широко используется. Несмотря на это, оно ограничено в том смысле, что 1) возбуждение должно быть четко определенной функцией и 2) дифференциальное уравнение должно иметь решение в замкнутой форме для данной функции возбуждения. Что, если возбуждение не следует четко определенной функции или у вас есть только случайные данные ввода-вывода? Это будет самый общий случай. Например, взгляните на график ниже, показывающий экспериментальные данные по напряжению и току для рассматриваемой нагрузки RL. Форма волны напряжения напоминает ступенчатую функцию, но определенно не идеальна. Мы могли бы просто аппроксимировать, что напряжение является ступенчатой ​​функцией, и работать с ней, но я думаю, что мы можем добиться большего.

Мы можем распространить описанный выше процесс на произвольные возбуждения, если вместо подгонки аналитического решения к данным мы подгоним смоделированный отклик системы на данные. Диаграмма ниже поясняет, что я имею в виду. Мы применяем одинаковую форму сигнала напряжения как к оборудованию, так и к нашему моделируемому оборудованию и наблюдаем за выходом каждого из них. Затем мы можем настраивать параметры системы (R и L) в моделировании до тех пор, пока моделируемый выходной сигнал не будет соответствовать измеренному отклику оборудования. Эту настройку можно даже автоматизировать с помощью процедуры оптимизации. И вишенка на торте — Python упрощает весь процесс.

Калибровка модели с помощью Python

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

Код

Как обычно, мы начнем с чтения всех необходимых модулей. NumPy, Pandas и Matplotlib довольно стандартны, но мы также импортируем несколько специфических функций из SciPy (назначение этих функций станет ясно позже). Пока мы этим занимаемся, мы установим некоторые настройки построения графиков по умолчанию в Matplotlib, чтобы сделать все наши графики согласованными и визуально привлекательными.

Экспериментальные данные хранятся в файле .csv со столбцами времени, напряжения и тока. Мы прочитаем данные в Pandas DataFrame и посмотрим на первые несколько строк, используя метод head.

Цифры никогда не бывают такими информативными, как реальная визуализация… Давайте нанесем на график данные о напряжении и токе, чтобы напомнить себе, с чем мы работаем.

Хорошо, вот тут-то и начинается самое интересное. Ранее мы говорили о том, что нам нужно подать одинаковое напряжение как на реальное, так и на смоделированное оборудование, чтобы мы могли сравнить их отклики. Чтобы смоделировать аппаратное обеспечение, мы собираемся создать класс RL_Load — если вы новичок в объектно-ориентированном программировании, класс можно рассматривать как план/шаблон для создания «объектов», которые в основном представляют собой просто контейнеры, которые хранить данные (называемые атрибутами) и иметь некоторую встроенную функциональность (называемую методами). Наш класс будет хранить параметры схемы как атрибуты и иметь несколько методов, которые позволят нам использовать их для моделирования.

Вкратце, но необходимо отметить, что моделирование обычно выполняется с использованием числового интегратора той или иной формы — эти числовые интеграторы требуют в качестве входных данных функцию следующего вида:

То есть им нужна функция, вычисляющая производную переменной состояния (в нашем случае текущей) по времени и текущему состоянию системы. Для нашей системы мы получили бы это, переписав уравнение цепи с точки зрения производной тока следующим образом:

Интуитивно приведенное выше выражение можно понять так: «Вы сообщаете мне текущее значение тока и какое напряжение приложено к цепи, а я могу сказать вам, как быстро меняется ток».

Хорошо, вкратце, давайте вернемся к коду. Глядя на определение класса ниже, вы можете видеть, что мы определяем несколько методов — мы рассмотрим каждый из них. Любой, кто знаком с объектно-ориентированным программированием на Python, узнает, что метод __init__ отвечает за инициализацию объекта. В нашем случае это просто сохранение параметров схемы. Далее метод di_dt делает именно так, как вы думаете; он вычисляет производную тока с учетом приложенного напряжения и текущего значения тока. Метод set_applied_voltage принимает функцию, которая возвращает приложенное к цепи напряжение как функцию времени, и сохраняет эту функцию как атрибут с именем V. Наконец, метод __call__ имеет особое значение в Python — его реализация позволяет нам вызывать наш объект, как если бы он был функцией. В нашем случае мы определяем метод таким образом, что вызов объекта RL_Load использует приложенную функцию напряжения для вычисления и возврата производной тока как функции тока и времени.

Чтобы сделать вещи более конкретными, в следующем фрагменте кода показано, как на самом деле используется класс. Сначала мы создаем объект RL_Load (со случайными значениями индуктивности и сопротивления) и назначаем его переменной rl. Затем мы создаем ступенчатую функцию для напряжения и «применяем» ее к rl, используя метод set_applied_voltage. Наконец, мы вызываем объект со входами время = 1 секунда и ток = 0 ампер и видим, что мгновенная скорость изменения тока в цепи в этот момент будет 1000 А/с.

Вы можете быть удивлены, но на данный момент мы в одной строке кода от запуска симуляции! Мы можем использовать функцию solve_ivp из SciPy для численного интегрирования нашей системы вперед во времени (обратите внимание, что ivp в solve_ivp означает «задача начального значения»).

Полную документацию о том, как использовать функцию solve_ivp, можно найти здесь, но ключевые входные данные: 1) функция, которая возвращает производную тока как функцию времени и тока (помните наше короткое отступление ранее?), 2) начальное и конечное время, между которыми мы хотим смоделировать, и 3) начальное значение тока. Глядя на приведенный ниже код, мы передаем нашу переменную rl непосредственно в интегратор (это можно сделать, потому что мы сделали объект вызываемым с помощью метода __call__). Результатом solve_ivp является объект, содержащий решение из моделирования (как а также несколько других показателей моделирования). Мы можем извлечь векторы решения времени и тока из этого объекта, а затем построить его, чтобы увидеть реакцию нашей созданной RL-нагрузки на входное напряжение единичного шага.

Мы собираемся сделать несколько графиков, подобных приведенному выше, поэтому вместо того, чтобы каждый раз переписывать код, мы создадим функцию, которая принимает объект RL_Load и экспериментальные данные (как Pandas DataFrame) и отображает смоделированные и измеренный текущий отклик на том же графике. Обратите внимание, что в вызове solve_ivp ниже время окончания симуляции основано на продолжительности экспериментальных данных.

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

Давайте соберем все до этого момента вместе. Приведенный ниже код создает объект RL_Load, применяет к нему экспериментальное напряжение, моделирует текущую реакцию, а затем отображает смоделированные и измеренные данные для сравнения.

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

Мы можем получить первоначальную оценку сопротивления, взглянув на отношение напряжения к току в установившихся условиях (т.е. ближе к концу измеренных данных). Кроме того, мы можем видеть, что кривой требуется около 0,035 секунды, чтобы достичь устойчивого состояния — зная, что для достижения устойчивого состояния системе первого порядка требуется около 4–5 постоянных времени (при условии, что входное напряжение на самом деле является ступенчатой ​​функцией). , мы можем оценить постоянную времени цепи.

Используя предполагаемое сопротивление и постоянную времени, мы можем опровергнуть оценку индуктивности, отметив, что постоянная времени цепи определяется отношением индуктивности цепи к сопротивлению (L/R). Если предыдущее обсуждение не было полностью ясным, не волнуйтесь, поскольку детали не так важны — ключевой момент, который я пытаюсь сделать здесь, заключается в том, что обычно мы можем быть более систематичными в нахождении начальной оценки для наших параметров. чем просто случайное угадывание.

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

На этом этапе мы могли просто продолжать угадывать и проверять различные значения параметров схемы, пока смоделированная кривая не соответствовала измеренным данным в достаточной степени. Однако это не очень хорошо обобщается. Например, что, если бы наша система имела 8 настраиваемых параметров?… или даже больше?! Не обязательно будет очевидно, как каждый параметр влияет на общее качество подгонки, не говоря уже о том, насколько утомительным будет процесс настройки. В качестве лучшей альтернативы мы будем использовать оптимизацию, чтобы найти «лучшую» оценку параметров схемы.

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

Метод __init__ берет экспериментальные данные, создает непрерывные функции интерполяции для напряжения и тока и устанавливает количество точек сетки, в которых оценивается ошибка между смоделированным и измеренным откликом. Метод __call__ делает вызываемым класс CostFunc (как мы обсуждали ранее) и отвечает за фактическую оценку того, насколько хороши оценки параметров. Он принимает вектор x (содержащий параметры схемы и начальное значение тока), создает и моделирует объект RL_Load с использованием входных параметров, а затем оценивает и возвращает среднеквадратичную ошибку между смоделированными и измеренными данными.

Мы включаем начальное значение тока в качестве настраиваемого входного параметра в оптимизацию, чтобы учесть тот факт, что наши измеренные данные тока могут быть не полностью обнулены в начале.

В следующем фрагменте показан конкретный пример использования класса CostFunc. Объект CostFunc создается с использованием экспериментальных данных и присваивается переменной cost. Затем мы вызываем объект CostFunc, передавая значения индуктивности, сопротивления и начального тока, и он возвращает скалярное значение, оценивающее, насколько «хороши» параметры.

На данный момент мы готовы к оптимизации. Существует несколько библиотек Python, поддерживающих оптимизацию, и широкий спектр алгоритмов, которые можно использовать, но мы будем использовать функцию minimize из SciPy с алгоритмом Нелдера-Мида. Полную документацию по функции minimize можно найти здесь, но ключевые входные данные включают функцию стоимости и начальное предположение для настраиваемых параметров. Результатом minimize является объект, содержащий оптимизированное решение и различные другие метаданные.

Обратите внимание, что выбор использования алгоритма Нелдера-Мида несколько произволен, но если вам интересно узнать больше, вы можете прочитать о каждом из доступных алгоритмов в minimize функции документация.

Этот следующий шаг не является абсолютно необходимым, но мы поместим функцию minimize в нашу собственную функцию optimization. Функция примет исходные оценки параметров схемы и экспериментальные данные, запустит оптимизацию, а затем распечатает и вернет оптимальные параметры схемы (в ожидании успешного завершения оптимизации). После определения функции мы запускаем оптимизацию с начальными оценками параметров, вычисленными ранее, и видим, что оптимизация успешно нашла набор лучших параметров!

Последний шаг — просто проверить, насколько хорошо он подходит. Приведенный ниже код создает новый объект RL_Load с использованием оптимизированных параметров схемы, а затем моделирует его для сравнения с экспериментальными данными. Эврика! Подгонка фантастическая, и у нас есть хорошая оценка для неизвестных параметров нашей схемы!

Заворачивать

Калибровка модели дорога моему сердцу, и я думаю, что многие инженеры и компании выиграют от ее более широкого использования. С относительно небольшим количеством кода мы смогли систематически оценивать неизвестные параметры RL-цепи, используя только данные о напряжении и токе.

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

Не стесняйтесь оставлять любые комментарии или вопросы, которые могут у вас возникнуть, или связаться со мной в Linkedin — я буду более чем счастлив прояснить любые моменты неопределенности. Наконец, я призываю вас поиграть с кодом самостоятельно (или даже использовать его в качестве начального шаблона для ваших собственных рабочих процессов) — код и соответствующие данные для этой статьи можно найти на моем Github.

Николас Хеменуэй

Больше контента на plainenglish.io. Подпишитесь на нашу бесплатную еженедельную рассылку новостей. Получите эксклюзивный доступ к возможностям написания и советам в нашем сообществе Discord.