Идентификация параметров математической моделиплоского движения сустава тема автореферата и диссертации по механике, 01.02.01 ВАК РФ

* #

Московский Государственный Университет имени М.В.Ломоносова

Механико-математический факультет Кафедра прикладной механики и управления

На правах рукописи УДК 531/533+576.72

Мартьянова Ольга Вениаминовна Идентификация параметров математической модели

плоского движения сустава (Специачьность -01.02.01-теоретическая механика)

АВТОРЕФЕРАТ

диссертациии на соискание учёной степени кандидата физико-математических наук

Москва-1994

Работа выполнена на кафедре прикладной механики и управления механико-математического факультета*Московского Государственного университета им. М. В Ломоносова.

Научный руководитель: доктор физико-математических

наук, профессор Е.А.Девянин

Официальные оппоненты: доктор физико-математических

наук, профессор Ю.Ф.Голубей кандидат физико-математических наук А.В.Ленский .

Ведущая организация: Институт проблем передачи

информации РАН

Защита диссертации состоится" IТ" ф е § р& ля!995г.

В |6 часов на заседании диссертационного Совета Д.053.05.01 при Московском Государственном университете по адресу: Москва 119899, Ленинские горы, Главное здание МГУ, корпус А, ауд.. 16-10

С диссертацией можно ознакомится в библиотеке механико-математического факультета МГУ.

Автореферат разослан "16 " 1995г.

Ученый секретарь диссертационного Совета: доктор физико-математических наук Л,,Е>. Тр&Ш&'Ь

I. Общая характеристика работы.

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

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

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

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

которой набор идентифицированных параметров определялся

однозначно. *

Методы исследований. В работе используются:

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

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

экспериментально-вычислительную установку. На

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

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

Аппробация. Результаты работы были доложены и обсуждены на семинарах кафедры прикладной механики и

управления механико-математического факультета МГУ, на семинаре "Роботы с элементами искусственного интеллекта" (под рук. академика РАН Д.Е.Охоцимского), на научной конференции "Биомеханика-94" (Москва, февраль 1994), на юбилейной научно-практической конференции ЦНИИПП-50, СПбНИИП-75 (Москва, май 1994), на второй всероссийской конференции по биомеханике (Нижний Новгород, ноябрь 1994)

Публикации. По теме диссертации автором опубликовано три работы, список которых приводится в конце автореферата.

Структура и объем диссертации. Диссертация состоит из введения, трех глав, приложения и списка литературы, .80 наименований. Общий объем диссертации 137 страниц.

II. Обзор содержания диссертации.

л

Во введении излагается краткая история данной задачи и приводятся основные результаты работы.

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

Активная мышца содержит сократительный элемент (СЕ), последовательный (БЕ) и параллельный (РЕ) "упругие" элементы.

Здесь СЕ - сократительный элемент, активное состояние которого регулируется нервно-мышечной системой, л$Е -последовательный пассивный упруго-вязкий элемент (сухожилие), РЕ - параллельный пассивный упруго-вязкий элемент.

Для описания элементов модели мышца рассматривается как макрообъект -. упруго-вязкое тело, обладающее двумя степенями свободы. Вводятся координаты ^ ,х, у как длины соответствующих элементов, ^ = х+у (см. рис.) fPE^fSE^fCE,

силы, развиваемые этими элементами. В общем случае они определяются как функции координат и t Kopouvii »н'мнпоп. Л сила, развиваемая сократительным элементом, зависит ещё и от управляющей переменной Е.

Ориентируясь на задачу идентификации, автором предложена следующая модель сустава:

Ja=Shi(a)^pi4i(a)+kd|i(a)-xi^vàcosa+MBH

©I 1=1 / . ^

Здесь:

а - суставный угол (см. рис. с tï>.?)

Х1 = х}-Х<н-относительная длина сократительного элемента i-ой мышцы

Xj -длина сократительного элемента i-ой мышцы

Xoi -длина сократительного элемента i-ой мышцы в

ненапряженном состоянии

hj -плечо i-ой мышцы в зависимости от а;

kpi ,ksi - коэффийиенты жесткости параллельного и

последовательного элементов i-ой мышцы;

Çi(a) = Çi(a)-I0i -относительная длина i-ой мышцы в

зависимости от а.

£,i(a) - общая длина i-ой мышцы в зависимости от а.

loi " общая длина i-ой мышцы в ненапряженном состоянии

J - момент иперппи суппнп;

Е1 -параметр,описывающий напряжение ьой мышцы ф|(Х|) -функция от Х| .описывающая максимальную

силу,которую может развивать 1-я мышца при длине сократительного элемента,равной

V] -параметр,характеризующий максимальную скорость

сокращения ьой мышцы V -вязкость.

Рассматривается модель с упрощенной геометрией сустава, (см. рис.), где срКх}) - квадратичная функция от Х^ Аь В|, С| -

коэффициенты квадратичного трехчлена. Тогда имеем:

И1(а)=Ь1со5а ^<0 Ь2 (а)= Ь2 соэа Ь2 > О ^(а^Цяла ^(а)=И25ша

вн

(

кД^ вта-Х; )=Е! + +С1) 1+- 0=1,2)

< i:, к,|—ЛЛЛЛ— AAAA SI

Kpi

VI

АЛАЛ-

KP2

V2

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

и=а

0 = _ У^яК^ Бта-х^^ (х)/х+ф;(х)]

г дй Ксоз2а к«х,Ь1+кц7х-,Ь')ч . уивта ^вн О „ =—=-—-+(~Ш-1—Т ) 81па+---+—а-

J

cosa

п да J

^ _ = дА = ^[к

' Еф^Х;)

Доказаны следующие утверждения.

Утверждение 1. Если (й,а,Х|,х2) положение равновесия, то

Dнi<0.

Утверждение 2: Если ¿п<0 , то положение равновесия

устойчиво. (

. _ с -(Ь,0п,+Ь,0Р9)

Утверждение 3: Если О < ——^-*—, то положение

п Г) Г)

равновесия неустойчиво.

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

Для идентификации параметров модели предлагается следующий набор опытов.

1 -ый опыт.

Обе мышцы расслаблены полностью. Е} = Е2 = О

Опыт состоит из серии статических измерений, т.е. .а=а=х=0

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

мвн(а)

Система уравнений при этом вырождается в алгебраическое выражение:

2-ой опыт.

Сгибатель напряжен, разгибатель полностью расслаблен. Угол ОС - фиксированная константа. Проводится серия подопытов для различных 0(, начиная с а = -20° и до а = +20° с шагом «2°. Во время каждого подопыта жестко фиксируется угол. Испытуемый начинает плавно напрягать сгибатель, пока не разовьет заданную силу. Фиксируются показания тензодатчика (т.е. развиваемый момент) и параллельно записываются миограммы обеих мышц. Миограммы интегрируются. Получаем

серию зависимостей момента от миограммы. Далее берется mio=0.5 и строится зависимость момента от (%.

3-ий опыт

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

Уриииснин дли »торою он 1»па принимают инд:

a=à=0, Ej =05 ,Е2 = 0, х, = О < М„„ =0

k^iicO-x^E^A.xf+B.x^C,)

4-ый опыт.

Свободное качание полностью расслабленной ноги.

е,=е2=о

Получаем следующие уравнения:

ч 1а=ЕЬ;(а)[кр;^(а)+к5}^(а)-х;)]-У(хсоза+Мвн

^(¿(a)-Xi)=0 (i=l,2)

снимаются показания (А и M(t).

5-ый и аналогичный ему 6-ой опыт.

Это опыт быстрой разгрузки. Фиксируется определенный угол Of 0 Испытуемый напрягает сгибатель. Для контроля

записывается параллельно миограммы обеих мышц и показания тснзодатчика. Резко отпускается фиксатор угла, и записывается зависимость 0C(t) .Причем, выбирается тот отрезок зависимости, где Е, = const (т.е. интегрированная миограмма представляет собой константу), а Е2 близка к нулю. Опыты, где не удалось

этого достичь, отбраковываются. Получаем уравнения:

-х^усссоБа+М

вн

(

Здесь необходимо производить отбор опытов, те отбраковать подопыты, где, например, миограмма антагониста недостаточно мала или где кривая момента недостаточно гладкая.

Данная последовательность опытов обладает следующими свойствами.

1) Все опыты являются вырождением общего случая системы (1), что существенно упрощает задачу идентификации.

2) Опыты упорядочены по усложнению вида системы (1), т.е. в последующих опытах можно использовать значения параметров, добытые из предыдущих опытов.

3) В данном наборе отражены основные характерные движения суставов, и, в то же время, эти движения являются простейшими.

4) Подобные опыты использовались в том или ином виде большинством исследователей, и, в этом смысле, этот набор опытов является классическим.

5) Опыты ориентированны на отслеживание параметра Е=соп51, т.к. вопросы управления "уровнем возбуждения" мышцы представляются слишком сложными для математического моделирования.

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

Первый этап метода идентификации состоит в. получениии теоретической зависимости. В случае динамических опытов теоретическая зависимость - решение дифференциальных уравнений с произвольно взятыми константами. В случае статических опытов теоретическая зависимость получена в виде уравнений. При решении дифференциальных уравнений применяется модифицированный метод Эйлера. Для того, чтобы сравнить полученную теоретическую зависимость с экспериментальной, необходимо вычислить экспериментальную зависимость именно в тех точках, в которых была получена теоретическая. Для этого применяется метол сглаживания функций, заданых поточечно, который можно назвать модификацией метода сплайнов третьего порядка. Второй этап является общим как для статических, так и для динамических опытов. Он состоит в следующем. Вычисляется сумма квадратов разностей между экспериментальной и теоретической зависимостями, т.е. численный аналог расстояния между функциями в метрике Ь2: Г = %)(х) - g2(x)¡ ёх. Будем рассматривать это расстояние как многопараметрическую функцию Г(с1,с2,...,ся), где С1,Сг,...,См- параметры, подлежащие идентификации в данном опыте. Задача идентификации параметров может быть интерпретирована как задача нахождения таких значений параметров, при которых функция принимает минимальное значение. Далее вычисляется

,, [д£ дТ аМ

градиент функции ^ =

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

величины смещения происходит следующим образом. Вначале устанавливается фиксированный коэффициент а=чао>0 и на

пересчёт градиента не производится, если при движении в том же направлении с тем же смещением значение функционала продолжает уменьшаться. Успешно пройдя некоторое число шагов, коэффициент а увеличивается вдвое и пересчитывается градиент, чтобы скорректировать направление смещения. Как только па некотором шаге значение функционала Г(с1,С2,...,См) оказывается больше, чем на предыдущем (это возможно уже вблизи локального минимума), то смещение п этом направлении не производится, а уменьшается вдвое, и пересчитывается градиент.

Когда значение коэффициента а становится меньше некоторого порогового значения (в программе пороговое

находимся вблизи локального минимума. Кроме того, в программе, реализующей этот метод, существует возможность вручную прервать процесс идентификации, наблюдая за значением функционала Г.

Кроме того, приближаясь к локальному минимуму, следует избегать ситуации, когда при численном подсчёте градиента величина с!с, оказывается настолько велика, что точки +с!с;,...,ск) и находятся по разные

стороны от локального минимума. В этом случае 1-ая компонента вычисленного таким образом градиента окажется противоположной по знаку ьтой компоненте истинного градиента. Поэтому. одновременно с уменьшением

первом шаге смещение равно

следующем шаге

значение равно

то можно сделать вывод, что мы

коэффициента а производится уменьшение и дна раза всех с^ 0 = 1,2,...,14) .

Преимуществами данного метода являются:

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

2.Высокая точность. Уменьшение величины смещения ах^ас)^ и величин с!с|, применяемых при численном подсчёте

градиента, обеспечивает достижение любой наперёд заданной точности вблизи локального минимума. Эта точность ограничена числом разрядов, выделяемых в программе для коэффициентов а, с1с{ и компонент градиента.

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

4.Универсальность. Алгоритм достаточно хорошо рлЬотлсг на широком классе задач: при различном числе идентифицируемых параметров, различных функционалах Г, различных требованиях на точность идентификации.

К недостаткам данного метода следует отнести то, что он не гарантирует нахождение глобального минимума в случае, когда существует несколько локальных минимумов, в достаточно большой окрестности которых градиент направлен в их сторону. Преодолеть этот недостаток можно либо варьируя начальную точку Г(Соьс025-">Соы)> либо теоретически доказав, что все

локальные минимумы такого рода находятся в ограниченной области, п II итой же о(шастн паходшен пюОллми.щ минимум. I) таком случае ошибку идентификации параметров С; можно оценить сверху размером данной области по ¡-ой координате.

.Третьим этапом является оценка точности идентификации параметров. После того, как найдена точка локального минимума (возможно, одного из нескольких вблизи глобального), исследуется его окрестность. Дли лого н шаре радиуса Г с центром в точке (С|,С1,...,См) строится £ -сеть и вычисляется значение функционата Г во всех узлах этой сети. Значение Г увеличивается до тех пор, пока указанный шар не будет заведомо включать точку глобального минимума (что контролируется с помощью значений функционала в узлах сети). Полученное Г и есть погрешность идентификации параметров.

Как правило, удаётся снизить погрешность для некоторых параметров С;, сузив шар до декартова произведения отрезков:

[С1 - ГЬС1 + п]х [С2 - г2,с2 + г2]х..^[ск -гк,сы + 1к].

В этом случае Г; даёт верхнюю оценку погрешности идентификации ¿-ого параметра.

Кроме того, во второй главе проанализировано поведение алгоритма идентификации для каждого из опытов: сходимость, оценка точности - радиуса Г шара, где находятся все локальные минимумы функционала Г.

В третьей главе производится анализ экспериментальных результатов и возникающих ошибок.

Ниже приводится таблица значений параметров, полученных в результате идентификации:

1.1 1.2 Ф) 2 3

J 0.558 0.559 0.732 0.596 0.780

КР1 192 192 193 196 286

Кв1 2910 2830 2570 2625

К£2 2600 2860 2520 2830

А1 . -2560 -2980 -3960 -3190 -2100

А2 -3420 -3100 -4380

В1 -2890 -3120 -3190 -3140

В2 -3130 -3120 -2810 -3280

С1 254 247 293 261

С2 245 245 399 233

VI 0.164 0.164 0.162 0.169 0.170

VI 1.92 1.86 2.91 2.33

Здесь приведены результаты пяти опытов с тремя различными испытуемыми.

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

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

сокращения от внутренних свойств мышц и сухожилий, что естественно было ожидать. <

Все возникающие ошибки можно разделить на несколько классов:

1. Ошибки измерения, связанные С тем, что все измерительные приборы, а именно, 'потенциометр, миографы и тензодатчики, имеют свою точность измерения.

2. Ошибки АЦП. Допустим, что сигнал, поступающий на АЦП в течение эксперимента, меняется в диапазоне [Р1,Р2]. АЦП преобразует его в целые числа из интервала [М1№|, т.е. при преобразовании возникает ошибка, равная (Р2-Р1)/(№-М1).

З.Ошибки калибровки. Во время калибровки происходят колебания системы, так, например, при небольшом изменении посадки испытуемого смещается положение равновесия для суставного угла, поэтому, данные, поступающие с АЦП усредняются ,что, естественно отличается от ситуации, если бы система была строго неподвижна. Таким образом возникают ошибки усреднения. Колебание напряжения в сети сказывается на том, что сигнал от датчиков так же колеблется, что оказывает своё влияние как при калибровке так и при последующих измерениях.

4.0шибки, эксперимента, специфические для различных опытов.

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

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

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

- в четвёртом, пятом и шестом опытах мы пренебрегаем трением в шарнирах механизмов.

5. Ошибки алгоритма.

- ошибки при решении дифференциальных уравнений

- ошибки аппроксимации функций, заданых поточечно.

- ошибки идентификации при нахождении локального минимума.

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

В приложении дано описание экспериментальной установки, приведены тексты программ и результаты опытов.

Основные результаты работы заключаются в следующем:

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

рассмотренной модели и модели, предложенной Д.Г.Алиевой, с точки зрения возможности идентификации параметров.

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

3.На экспериментальной установке проведены серии опытов с несколькими испытуемыми, для каждого из которых определены конкретные параметры выбранной математической модели.

4. Проведен сравнительный анализ полученных данных для различных испытуемых. Проверена возможность модели с идентифицированными параметрами предсказывать характер экспериментальных зависимостей.

5.Проведен аналнз возникающих ошибок измерений и погрешностей вычислений.

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

Публикации автора по теме диссертации

1.Мартьянова О.В. Об идентификации параметров уравнения движения сустава под действием пары мышц-антагонистов. - В сб.:"Вопросы проектирования информационных и кибернетических систем". Уфа 1992.

2.Девянин Е.А., Мартьянова О.В. Движение сустава, обеспечиваемое парой мышц-антагонистов,- Тез. докл. на юбилейной научно-практической коференции ЦНИИПП-50,СПбНИИП-75. Москва, 1994 г.

3.Мартьянова О.В. Исследование плоского движения сустава под действием мышц-антагонистов. 'Ici. Докл.на второй всероссийской конференции по биомеханике. Нижний Новгород, 1994 г.