Санкт-Петербургский Государственный Электротехнический Университет «ЛЭТИ» Кафедра БТС Моделирование в биотехнических системах Анна Юрьевна Виллевальде 2011
Моделирование – исследование объектов познания на их моделях. Моделирование – неотъемлемая часть всякой целенаправленной деятельности. Модели: • Цели • Планы • Представления 2
Мир моделей субъекта Модели: • Врожденные • Приобретенные Реальный мир 3
Модель – заменитель оригинала • Подобный оригиналу в конечном числе отношений • Более удобный в использовании, чем оригинал Предпосылки использования моделирования - формирование теории систем (методологии изучения сложных систем); - мощное развитие вычислительной техники, ее доступность => накоплен мощный научно-технический потенциал в виде программных продуктов, в том числе касающийся реализации системных принципов. 4
Модель (лат. modulus – мера) – это объект (материальный, представление об объекте, математическая зависимость, программа для компьютера), находящийся в отношении подобия к моделируемому объекту. – это объект-заместитель объекта-оригинала, инструмент для познания, который исследователь ставит между собой и объектом и с помощью которого изучает некоторые свойства оригинала. В качестве модели выступает другой материальный или мысленно представляемый объект, замещающий в процессе исследования объект-оригинал. – это некоторая вспомогательная искусственная или естественная система, находящаяся в некотором объективном соответствии с познаваемым объектом, способная замещать его в определенных отношениях, дающая при её исследовании, в конечном счете, информацию о самом моделируемом объекте. 5
Метасистема моделирования Культура (среда моделирования) Модель Субъект Объект 6
Модель – Объект Модель – конечное отображение оригинала • Модель – не сам объект; • Ресурсы моделирования ограничены; • Происходит отбор свойств. Модель подобна оригиналу лишь в конечном числе отношений. Модель – упрощенное отображение оригинала Упрощения необходимы и допустимы: • Упрощение – сильное средство выявления главной закономерности; • Снижение затрат на моделирование. 7
Модель – Объект Модель – приближенное отображение оригинала • Модели подобны оригиналу в разной степени приближенности; • Возможность увеличивать степень приближенности к оригиналу. Модель – целевое отображение оригинала • Для модели отбираются только те свойства, которые отвечают цели моделирования; • Могут быть разные и даже противоречивые модели для одного объекта. 8
Модель – Объект Истинное и ложное в модели I II III Отсутствует в модели: – Неизвестно; – Известно, но несущественно; – Известно и существенно, но не включено по ошибке. Истинное – Известно точно (достоверно); – Известно с некоторой неопределенностью; – Достоверно при выполнении условий. Безусловно ложное – Существенно не влияет на результат; – Порождает ошибки моделирования. Объект I Модель II III Модель имеет как безусловно истинное, так и условно истинное и безусловно ложное содержание. ИСТИННОСТЬ • Сама модель не может быть истинной или ложной; • Степень истинности выявляется при соотнесении с реальностью; • Важны условия соотнесения! 9
Модель – Субъект Адекватность модели Модель, с помощью которой успешно достигается поставленная цель, называется адекватной • Познавательные модели (пассивная деятельность); М форма организации и представления знаний, средство соединения новых знаний с уже имеющимися. • Прагматические модели (активная деятельность): средство управления, средство организации практических действий, способ представления образцово правильных действий или их результата. • Для познавательных моделей: истинность = адекватность. • Для прагматических моделей: адекватной может быть ложная модель. • Могут быть разные адекватные модели. О 10
Ингерентность моделей Модель – Культура – достаточная степень согласованности создаваемой модели со средой, в которой ей предстоит функционировать, чтобы она входила в эту среду не как чужеродный элемент, а как естественная составная часть (в модели должны быть предусмотрены «стыковочные узлы» со средой, интерфейсы). Другой аспект ингерентности модели: в самой среде должны быть созданы предпосылки, обеспечивающие функционирование модели. Т. е. не только модель должна приспосабливаться к среде, но и среду необходимо приспосабливать к модели. Модель может выполнять свои функции только в зависимости от обстоятельств. 11
Динамика моделей Жизненный цикл моделей: Жизненный цикл моделей • Возникают, изменяются, взаимодействуют, уступают место другим моделям; • Могут жить дольше субъекта. Трудности моделирования • • Модель в разной культуре; Противоречивость требований; Цели могут меняться; Построение модели – всегда сложно. 12
Моделирование Материальное (предметное, реальное) и идеальное (информационное, абстрактное) Материальная модель – некоторый материальный объект. Реальные модели создаются из реальных объектов. Физические (вещественные, натурные) и аналоговые (энергетические) модели 13
Материальные модели Физические модели воспроизводят структуру исследуемого объекта и взаимоотношения его частей: уменьшенная или увеличенная копия объекта, допускающая его исследование с последующим переносом свойств модели на реальный объект с помощью теории подобия. Напр. , макет здания, самолета, модели молекул веществ. Аналоговые модели не являются копией структуры исследуемого объекта, могут иметь иную физическую природу; формально воспроизводят свойства, связи и другие характеристики объекта, что позволяет экспериментировать с моделью и исследовать оригинал. Напр. , груз на пружинке – аналоговая модель многих колебательных процессов, включая биологические; эффективность биопрепаратов и других воздействий на БО исследуется на биологических моделях, на которых воспроизводится та или иная патология или реакция на воздействие: моделирование гриппа на белых мышах, лучевой болезни – на свиньях. 14
Материальные модели Свойства материальной модели подобны свойствам объекта за счет: • Прямого подобия: достигается за счет результатов физического взаимодействия или цепочки взаимодействий модели и оригинала. Напр. , фотография, макет. • Косвенного подобия: достигается за счет совпадения свойств моделей (модели-аналоги). Напр. , физическая аналогия, историческая аналогия. • Условного подобия: достигается за счет договоренности, условленности. Напр. , деньги, буквы, слова, знаки. 15
Моделирование Идеальное (информационное) моделирование основано не на материальном сходстве, а на аналогии идеальной, информационной, мыслимой. Интуитивное (вербальное, мыслительное) и знаковое (языковое) моделирование 16
Идеальные модели Интуитивные модели основаны на интуитивном представлении об ОИ, не поддающемся формализации или не нуждающемся в ней. Носят описательный характер. Напр. , жизненный опыт каждого человека может рассматриваться как интуитивная модель окружающего мира. Знаковые модели основаны на знаковых представлениях и преобразованиях какого-либо вида: схемах, графиках, чертежах, формулах, листингах программ, а также включают совокупность законов, по которым можно оперировать с выбранными знаковыми отображениями и их элементами. Важнейшим видом знакового моделирования является математическое моделирование. Иерархия языков: Иерархия языков - Разговорный - Профессиональный … - Математический Расплывчатость Точность 17
Математическое моделирование – это метод исследования процессов или систем путем построения их математических моделей и исследования этих моделей. – процесс установления соответствия данному реальному объекту некоторого математического объекта, называемого математической моделью, и исследование этой модели, позволяющее получать характеристики рассматриваемого реального объекта. Вид математической модели зависит как от природы реального объекта, так и от задач исследования объекта, а также требуемой достоверности и точности решения этих задач. 18
Математическая (аналитическая) модель – это система математических соотношений, описывающих изучаемый процесс или явление. – это «эквивалент» объекта, отражающий в математической форме важнейшие его свойства: законы, которым он подчиняется, связи, присущие составляющим его частям и т. д. – это совокупность математических соотношений, уравнений, неравенств и т. п. , описывающих основные закономерности, присущие изучаемому процессу, объекту или системе. – абстрактное математическое представление процесса, устройства или теоретической идеи; использует набор переменных, чтобы представлять входы, выходы и внутренние состояния, а также множества уравнений и неравенств для описания их взаимодействия. – уравнение, выражающее идею. 19
Имитационная модель – знаковая модель, имитирующая функционирование системы. В отличие от классических аналитических моделей, имитационные модели не обязательно описываются с помощью определенного математического аппарата. В процессе построения модели осуществляется переход от реальной системы к определенной логической схеме. При этом определяются компоненты или элементы системы, параметры, переменные, связи между элементами, функциональные зависимости. Последние могут быть представлены в виде самых разнообразных аналитических выражений. В конечном итоге модель компонуется в виде программы для компьютера. Описание модели на языке программирования объединяет в единое целое все элементы системы, «оживляет» их и дает возможность ставить машинные эксперименты с моделями, имитируя поведение систем. 20
Классификация моделей • • • статические – динамические; дискретные – непрерывные; детерминированные – стохастические; сосредоточенные – распределенные; линейные – нелинейные; стационарные – нестационарные 21
Классификация моделей • Модель называется статической, если фактор времени не учитывается или его значение зафиксировано (рассматривается временной срез). Модель называется динамической, если ее функция зависит от процесса (его предыстории), где включено время. В таких моделях может учитываться влияние памяти и инерционности процессов. Предметом изучения при моделировании является изменение рассматриваемого объекта во времени. • Дискретность / непрерывность зависит от того, дискретными или непрерывными являются множества входных и выходных величин модели и значений фактора времени. Если хотя бы одно из этих множеств дискретно, то модель относится к классу дискретных. Дискретные величины – это конечное множество величин, которые принимают «оторванные» друг от друга значения, допускающие нумерацию в числах натурального ряда. Непрерывная величина принимает все значения из некоторого заданного интервала. Напр. , моделями непрерывной величины являются отрезок, луч, прямая; дискретной моделью по входу и выходу является модель реакции мышечного волокна на импульсное раздражение, волокно может находиться в одном из трех состояний: покоя, возбуждения или рефрактерности (когда оно не чувствительно к 22 раздражению).
Классификация моделей Детерминированные модели – такие модели, для которых задание входа однозначно определяет выход модели. Стохастические модели – для которых задание входа неоднозначно определяет выход элементарные вероятностные модели в виде законов распределения случайных величин. Модели с сосредоточенными и распределенными параметрами относятся к дифференциальным уравнениям. Если модель описывает зависимость одного или нескольких неизвестных от одного аргумента, то такая модель – с сосредоточенными параметрами. Модель с распределенными параметрами обычно описывает зависимость неизвестной одновременно от нескольких аргументов, напр. , временной и пространственных координат. Напр. , зависимость температуры потока крови в сосуде от времени и расстояния до любой точки потока от начала движения. Обычно такие д. у. представляются в виде частных производных. Линейные / нелинейные модели определяются линейностью или нелинейностью связи входных и выходных параметров. В биологических системах преобладают нелинейные связи, что существенно осложняет их моделирование. Определение стационарности модели относится к моделированию случайных процессов. В стационарных моделях (в отличие от нестационарных) характеристики 23 не меняются во времени (спектральные характеристики,
Классификация моделей по глубине описания объекта (исходя из глубины отражения структуры и свойств описываемой системы) 1. Функциональные модели – описывают объект как черный ящик, т. е. о механизме функционирования которого ничего не известно, имеются лишь входные и выходные переменные и параметры. Такая модель не дает прямой информации о механизме функционирования системы, однако она позволяет управлять процессом, оптимизировать его, исследовать влияние комплекса факторов и т. д. Напр. , оценка и оптимизация влияния воздействия одного или нескольких лечебных препаратов на организм. 2. Структурно-функциональные модели – в них учитывается функционирование структурно-обособленных подсистем, при этом могут моделироваться и взаимосвязи между ними, как горизонтальные (на одном уровне организации), так и вертикальные, в соответствии с иерархической структурой системы. Напр. , в процессе исследования новых лечебных препаратов строится комплекс функциональных моделей, описывающих зависимости выходных характеристик различных подсистем организма, определенным образом связанных между собой. 24
Классификация моделей по глубине описания объекта (исходя из глубины отражения структуры и свойств описываемой системы) 3. Модели, учитывающие механизмы функционирования исследуемой системы – в таких моделях используются известные или гипотетические сведения о структурных и причинноследственных связях между отдельными элементами системы, в результате становится возможным активно вмешиваться во внутренние механизмы исследуемых явлений в процессе исследования модели. Напр. , математическая модель развития бактериального или вирусного инфекционного заболевания, лечение этих заболеваний препаратами с различными механизмами действия. 25
Общие алгоритмы моделирования Процесс исследования БО на основе методов математического моделирования начало Постановка задачи исследования Идентификация ОИ как системы Построение модели Исследование модели Принятие решения конец 26
Задачи моделирования 1. Компактная и формализованная запись представлений о наблюдаемых явлениях и процессах, характеризующих функционирование исследуемой системы В результате появляется возможность, базируясь на небольшом числе исходных постулатов, описать и исследовать широкий круг явлений. При синтезе моделей осуществляется перевод описываемого процесса с языка соответствующей науки на язык математического моделирования. В итоге эти процессы, которые на разных языках выглядят по-разному, приводятся к единой форме, что позволяет вскрывать общие структурнофункциональные характеристики систем и более глубоко проникать в сущность процессов и явлений. Т. о. , становится возможным заимствовать достижения в изучении аналогичных систем из тех наук, где достигнут больший прогресс в данном направлении. Напр. , динамические процессы ряда систем в физике, электротехнике, биологии могут исследоваться с помощью одинаковых по структуре моделей – д. у. 2 го порядка, описывающих перемещение груза, прикрепленного к пружине, изменение электрического напряжения на выходе звена, состоящего из емкости и сопротивления, а также динамику выходных характеристик БО при дозированном внешнем воздействии. 27
Задачи моделирования 2. Проверка гипотезы о функционировании системы В биологии и медицине нередко бывает сложно проверить гипотезу о функционировании системы в эксперименте на живом объекте. Моделирование позволяет на основе результатов эксперимента, проведенного с моделью, подтвердить, скорректировать или опровергнуть выдвинутую гипотезу. Напр. , эффект парадоксального действия лекарственных средств нередко наблюдается врачами разнообразные гипотезы. 3. Выдвижение новых гипотез о функционировании системы Напр. , при моделировании лечения инфекционного заболевания была выдвинута гипотеза о влиянии концентрации возбудителя в организме на эффективность борьбы с заболеванием; с помощью лечебного препарата концентрация возбудителя регулируется таким образом, чтобы, с одной стороны, стимулировать защитноприспособительные силы организма, с другой – не нанести существенный вред организму (оптимальная концентрация). 4. Исследование поведения систем в экстремальных условиях – моделирование поведения живого объекта в экстремальных условиях. Напр. , с помощью моделирования можно оценить действие различных лекарственных средств в случае отравления или радиационного поражения; можно моделировать действия оператора в экстремальных условиях (основание для 28 профессионального отбора).
Задачи моделирования 5. Управление исследуемым объектом, оптимизация его структуры, функции Напр. , модель, построенная на основе теории массового обслуживания, может быть использована для оптимального функционирования поликлиник или использования медицинской техники. 6. Прогнозирование процесса Напр. , оценка и прогнозирование развития эпидемии (ВИЧ, гепатит, туберкулез). 7. Систематизация и сопоставление исследуемых объектов Математическое моделирование позволяет проводить сопоставление исследуемых объектов по совокупности параметров, характеризующих их структуру и функции, в результате становится возможным существенно повысить эффективность как сравнительной оценки объектов, так и эффективность их систематизации. Напр. , оценка эффективности иммуностимуляторов на фоне терапии антибиотиками (сложная, многофакторная). 29
Идентификация объекта исследования как системы Классические науки, особенно биологические, традиционно оперировали простыми одномерными зависимостями типа «доза – эффект» с целью определения причинно-следственных связей. Однако исследование целостного объекта требует более сложных функциональных зависимостей, охватывающих одновременно различные стороны его функционирования. Поэтому при постановке задачи математического моделирования необходимо рассматривать любой объект, даже кажущийся простым, как сложную систему. Система – это совокупность элементов и отношений между ними, образующих единое целое. Систему можно представить как совокупность двух множеств, X -входов и Y-выходов, и отношений между ними, формализующих связь между входом и выходом. Элемент системы – объект, входящий в систему в качестве относительно самостоятельной единицы, неделимой относительно заданного уровня моделирования. 30
Идентификация объекта исследования как системы Общие свойства систем: целостность, структурированность, целенаправленность. Целостность (единство) означает, что система отделена от внешней среды и взаимодействует с ней через входы и выходы. Структурированность: система включает в себя множество подсистем, определенным образом связанных и взаимодействующих между собой. Для биологической системы характерна многоуровневая иерархическая организация (система бесконечна вверх и вниз). Это означает, что элементы исследуемой подсистемы могут рассматриваться как подсистемы более низкого уровня, и наоборот, подсистемы низкого уровня превращаются в элементы системы более высокого уровня. Существуют также горизонтальные связи между подсистемами (на одном уровне). 31
Идентификация объекта исследования как системы 1. Выявление причинно-следственных отношений (связей) в системе и сбор информации о механизмах явлений. 2. Определение системных параметров и переменных: входных, выходных, управляемых, неуправляемых, зависящих от времени. 3. Определение степени сложности организации системы. Параметр – это переменная величина, которой присваивается постоянное значение в рамках указанного применения (на период функционирования модели). 32
Схематическое представление системы параметров и переменных объекта исследования x 1 x 2 xk … u 1, u 2, …, un ; a 1, a 2 , …, al … y 1 y 2 ym x 1, x 2, …, xk – входные параметры, переменные, вход системы (не только внешние факторы, но могут быть и внутренние; то, что влияет на выход); u 1, u 2, …, un – переменные, характеризующие свойства системы и изменяющиеся в процессе исследования (моделирования) системы; a 1, a 2, …, al – параметры, характеризующие свойства системы, постоянные в процессе моделирования; y 1, y 2, …, ym – выход системы (иногда выходом может быть u, напр. , критические значения u). 33
Основные функциональные связи, исследующиеся в процессе моделирования 1. U = F(X, A, t) ui = fi(x 1, …, xk; a 1, …, al; t) 2. Y = F(X, U, A, t) yi = fi(x 1, …, xk; u 1, …, un; a 1, …, al; t) 3. чаще всего используется Y = F(X, A, t) yi = fi(x 1, …, xk; a 1, …, al; t) в наиболее простом варианте yi = fi(x 1, …, xk; t) и yi = fi(x 1, …, xk) – не зависит от времени Задача: Лечение гипертонической болезни кордипином и эналаприлом 34
Определение степени сложности и организации системы может играть важную роль в выборе класса модели; Сложность системы является качественной характеристикой, для которой не существует единого формального метода оценки. Показатель сложности должен характеризовать структуру и поведение системы, а также трудности, связанные с ее изучением, созданием и использованием. В качестве наиболее примитивной количественной меры сложности системы следует рассматривать её разнообразие, т. е. количество возможных состояний: Hm = log n, n – количество состояний, которые принимает система. Если выходные параметры системы непрерывны, необходимо сначала осуществить их дискретизацию с шагом Δy, тогда n=(ymax–ymin)/Δy, Δy выбирается, исходя из требования точности оценки и возможностей прибора (метода). Организация системы: R = 1 – (H / Hm), H – оценка энтропии системы: Pi – вероятности состояний системы. Сложность – максимально возможная энтропия системы. Когда вероятности состояний системы равны между собой, H = Hm. Когда H = Hm, R = 0. 35
Определение степени сложности и организации системы Классификация систем по степени сложности и организации По степени сложности: 0 ≤ Hm < 3, n = 20 – простые системы; 3 ≤ Hm < 6, n = 400 – сложные системы; Hm ≤ 6 – очень сложные системы. По организации: детерминированные 0, 3 < R ≤ 1 квазидетерминированные 0, 1 < R ≤ 0, 3 вероятностные 0 ≤ R ≤ 0, 1 Напр. , детерминированные системы осуществляют управление внутренними органами: подсистемы, управляющие функционированием сердечнососудистой системы, реакция зрачка на свет. Вероятностные системы определяют взаимодействие анализаторов человека и поведенческих реакций: реакция человека на один и тот же стимул может быть различной, она зависит от многих факторов, которыми исследователь не может управлять. Квазидетерминированные системы также могут осуществлять управление внутренними органами, при условии, если ведущую роль играет ЦНС; некоторые патологические процессы. Задача: Исследование реакции БО на пороговое электрическое раздражение 36
Оценка организации динамических систем Динамическая система: если каждый раз значение на выходе системы, при одном и том же входном значении, разное, т. е. зависит от того, в какой последовательности подавались входные значения. Если система последовательно переходит из одного состояния в другое в пределах заданного множества состояний, то оценку ее организации необходимо проводить поэтапно, причем каждый этап соответствует определенному уровню оценки: На первом этапе (уровне) оценивается вероятность появления каждого единичного состояния, независимо от динамики других состояний (статическая задача). Начиная со второго уровня оцениваются вероятности появления комбинаций состояний по 2 (второй уровень), по 3 (третий уровень) и так до объединения в одну группу всех k состояний. Количество состояний n на каждом i-ом уровне определяется как количество сочетаний из k по l: k – количество состояний, l – количество состояний в объединяемой комбинации. Сложность, энтропия и организация на i-ом уровне оценки: В конечном итоге строится зависимость R от уровня оценки. Задача: оценить организацию процесса, близкого к синусоидальному и представляющего 37 последовательную смену 4 х состояний.
Выбор класса модели в значительной степени зависит от уровня организации и сложности ОИ; относится к классификации по признакам: детерминированные – стохастические и линейные – нелинейные модели. Hm 12 6 33 32 31 3 23 22 21 13 12 11 0, 1 0, 3 1 R 11 – простые детерминированные линейные модели; 12 – квазидетерминированные простые линейные модели; 13 – простые вероятностные модели; 21 – сложные детерминированные модели, линейные и нелинейные; 22 – сложные нелинейные квазидетерминированные модели; 23 – сложные вероятностные модели; 31, 32, 33 – очень сложные детерминированные, квазидетерминированные, вероятностные модели. 38
Выбор класса модели С помощью построенной сетки в определенной степени можно систематизировать выбор математического аппарата моделирования: Для наиболее распространенных классов простых и сложных детерминированных моделей (11 и 21) адекватен следующий математический аппарат: - Экспериментально-статистические модели, основанные на алгебраических уравнениях. Применяются для аппроксимации результатов эксперимента. - Методы дифференциальных и интегральных уравнений. Предусматривают построение математических моделей динамических систем. Позволяют достичь более глубокой содержательности, чем при использовании алгебраических уравнений. - Методы теории управления и оптимизации. Позволяют анализировать биосистемы, разделяя их на объекты управления и управляющие подсистемы, и решать прикладные задачи управления биосистемами, оптимизируя их функции. - Теория конечных автоматов. Адекватна для моделирования динамики четко выраженных детерминированных систем с дискретными состояниями. - Булева алгебра. Применяется для статических моделей (закономерности переходов не зависят от динамики). 39
Выбор класса модели Для простых и сложных квазидетерминированных моделей (12, 22): Нелинейная алгебра и дифференциальные уравнения, а также уравнения, коэффициенты которых подчиняются некоторым законам распределения. Для простых и сложных вероятностных моделей (13, 23): Теория вероятностей, теория случайных процессов, алгебраические и дифференциальные уравнения для описания вероятностей и законов распределения, метод Монте-Карло с использованием генераторов случайных чисел, Марковские цепи, теория распознавания образов. Однако в реальных условиях при моделировании систем часто необходимо учесть многие свойства системы, которые одним математическим аппаратом не охватить, в то же время совместное использование нескольких математических аппаратов часто приводит к неразрешимым задачам, например, использование дискретного и непрерывного аппарата вместе. Для разрешения таких противоречий используются методы имитационного моделирования, где могут одновременно применяться различные математические аппараты. 40
Расчет коэффициентов модели Расчет коэффициентов основан на сопоставлении расчетных и экспериментальных данных. Расчетные данные получают с помощью модели, а экспериментальные – путем проведения экспериментов с ОИ, из литературных данных или исходя из гипотез. Расчет коэффициентов для статических моделей Линейные регрессионные модели Задача регрессионного анализа: Зная множество значений на входах и выходах системы, построить модель, т. е. определить функцию системы, по которой ее вход преобразуется в выход. Линейная одномерная регрессионная модель: черный ящик, имеющий один вход и один выход; зависимость между входом и выходом линейная или почти линейная. 1) Исследователь вносит гипотезу о структуре ящика Рассматривая полученные экспериментально данные, предположим, что они подчиняются линейной гипотезе: Y = A 1 X + A 0 41
Расчет коэффициентов для статических моделей Линейные регрессионные модели 2) Определение неизвестных коэффициентов A 0 и A 1 модели Метод наименьших квадратов Для каждой из n снятых экспериментально точек вычислим ошибку Ei между экспериментальным значением Yi. Эксп. и теоретическим значением Yi. Теор. , лежащим на гипотетической прямой A 1 X + A 0: Ei = (Yi. Эксп. – Yi. Теор. ), i = 1, …, n; Ei = Yi – A 0 – A 1 · Xi, i = 1, …, n. Ошибки Ei для всех n точек следует сложить. Чтобы положительные ошибки не компенсировали в сумме отрицательные, каждую из ошибок возводят в квадрат и складывают их значения в суммарную ошибку F одного знака: Ei 2 = (Yi – A 0 – A 1 · Xi)2, i = 1, …, n Цель: минимизация суммарной ошибки F за счет подбора коэффициентов A 0, A 1. Необходимо найти такие коэффициенты A 0, A 1 линейной функции Y = A 1 X + A 0, чтобы ее график проходил как можно ближе одновременно ко всем экспериментальным точкам. 42
2) Метод наименьших квадратов Чтобы минимизировать суммарную ошибку, найдем частные производные функции F по каждой переменной и приравняем их к нулю (условие экстремума): После раскрытия скобок получим систему из двух линейных уравнений: Решение имеет вид: Для нахождения A 0 и A 1 методом Крамера представим систему в матричной форме: 43
Расчет коэффициентов для статических моделей Линейные регрессионные модели 3) Проверка Чтобы определить, принимается гипотеза или нет, нужно: 1. рассчитать ошибку между точками заданной экспериментальной и полученной теоретической зависимости Ei и суммарную ошибку F 2. найти значение σ где F – суммарная ошибка, n – общее число экспериментальных точек. Если в полосу, ограниченную линиями YТеор. -S и YТеор. +S, попадает 68. 26% и более экспериментальных точек, то выдвинутая гипотеза принимается. В противном случае выбирают более сложную гипотезу или проверяют исходные данные. Если требуется большая уверенность в результате, то используют дополнительное условие: в полосу, ограниченную линиями YТеор. -2 S и YТеор. +2 S, должны попасть 44 95. 44% и более экспериментальных точек.
Расчет коэффициентов для статических моделей Линейные регрессионные модели 3) S = σ/sin(β) = σ/sin(90° – arctg(A 1)) = σ/cos(arctg(A 1)) Условие принятия гипотезы выведено из нормального закона распределения случайных ошибок, P – вероятность распределения нормальной ошибки. 45
Расчет коэффициентов для статических моделей Линейные множественные модели Функциональная структура системы имеет линейную зависимость, но количество входных сигналов, одновременно действующих на объект, равно m. Анализируется аналогично одномерной модели. Нелинейные регрессионные модели Многие нелинейные модели при помощи подстановок и переобозначений приводятся к линейной множественной модели. Полиномиальная множественная регрессионная модель Если система имеет 2 входа, а зависимость выхода от входов напоминает квадратичную, то целесообразно выбрать такую гипотезу: Y = A 0 + A 1 · X 1 + A 2 · X 2 + A 3 · X 1 · X 2 + A 4 · X 1 + A 5 · X 2 Обозначим: Z 1 = X 1 · X 2; Z 2 = X 1 · X 1; Z 3 = X 2 · X 2 Y = A 0 + A 1 · X 1 + A 2 · X 2 + A 3 · Z 1 + A 4 · Z 2 + A 5 · Z 3 46
Нелинейные регрессионные модели Мультипликативная регрессионная модель Y = A 0 · X 1 A 1 · X 2 A 2 · … · Xm. Am Прологарифмируем левую и правую части данного уравнения: ln(Y) = ln(A 0) + A 1 · ln(X 1) + A 2 · ln(X 2) + … + Am · ln(Xm) Обозначим: W = ln(Y), B 0 = ln(A 0), Z 1 = ln(X 1), Z 2 = ln(X 2), …, Zm = ln(Xm) Получим: W = B 0 + A 1 · Z 1 + A 2 · Z 2 + … + Am · Zm Обратная регрессионная модель Y = k/(A 0 + A 1 X 1 + … + Am. Xm) Заменим: W = 1/Y, ai = Ai/k Получим: W = a 0 + a 1 · X 1 + … + am · Xm Экспоненциальная модель Y = e. B 0 + B 1 X 1 + B 2 X 2 + … + Bm. Xm Прологарифмируем левую и правую части уравнения: ln(Y) = B 0 + B 1 · X 1 + B 2 · X 2 + … + Bm · Xm Выполним замену W = ln(Y) и получим: W = B 0 + B 1 · X 1 + B 2 · X 2 + … + Bm · Xm 47
Примеры применения нелинейных регрессионных моделей Степенная функция Y = A · XB – связывает параметры организмов животных различных видов Y с характеристиками тела X. Напр. , Y – частота обменных процессов, пульса, дыхания, чувствительность к вредным воздействиям; X – масса тела, вес мозга, площадь поверхности тела. В результате логарифмирования уравнение сводится к линейной модели. Существуют специальные атласы с графическим представлением таких функций для различных видов животных. Такие модели могут использоваться для переноса данных, полученных в экспериментах с животными, на человека. 48
Примеры применения нелинейных регрессионных моделей В качестве объекта исследования рассмотрим гипсовый слепок нижней челюсти человека. Проекция оттиска на плоскость YOX позволяет получить набор точек, задающих внешний (или внутренний) профиль челюсти. Полученную выборку координат точек {xi, yi} можно рассматривать как совокупность значений входа x и выхода y. Выдвинута гипотеза о том, что внешний профиль челюсти можно описать с помощью уравнения параболы: Следует найти такие значения коэффициентов b 0, b 1, b 11, которые обеспечат минимум функционала: Минимум функционала достигается при выполнении требований: Откуда: 49
Примеры применения нелинейных регрессионных моделей Традиционная графическая форма представления тональной пороговой аудиограммы является отображением значений порогов слышимости звука на отдельных частотах, измеренных по воздушному и костному проведениям. По оси абсцисс в логарифмическом масштабе откладывают значения частот fi = 125, 250, . . . , 16000 Гц, на которых проводят исследования, по оси ординат – потери слуха Yi = -15, …, 110 д. Б; Y = 5 д. Б. Нулевой линией отмечены средние потери слуха, измеренные у хорошо слышащих людей. При наличии патологий органов слуха на аудиограмме фиксируется увеличение потерь слуха на отдельных частотах. Точка П 1 характеризует потери слуха по костному проведению на частоте 250 Гц; точка П 2 – потери слуха по воздушному проведению на частоте 1000 Гц. 50
Примеры применения нелинейных регрессионных моделей Известно, что аудиограммы одного и того же пациента, снятые с некоторым временным сдвигом в течение одного дня (или нескольких дней), будут различаться. Учитывая эти объективные наблюдения, рассмотрим возможность построения формализованного описания кривых потерь слуха с целью создания модели порогов слуха. Вполне очевидно, что в качестве зависимой переменной y надо рассматривать потери слуха. Т. к. речь идет об аппроксимации плоских кривых, то число независимых переменных не должно превышать 1. Единственной независимой переменной x является частота испытательного сигнала. Так как ее значения всегда постоянны, то для упрощения модели в качестве x будем использовать порядковый номер частоты испытательного сигнала: 1, 2, …, 11. Для построения уравнения модели используем полиномиальное уравнение вида: m – порядок полинома. Для оценки коэффициентов уравнения используем выборку из нескольких аудиометрических кривых, полученных по воздушному проведению при исследовании слуха одного пациента: Каждый эксперимент представлен одной строкой, в последнем столбце приведены оценки дисперсии потерь слуха, в последней строке – средние значения потерь слуха на каждой частоте (Ср. ПС). 51
Примеры применения нелинейных регрессионных моделей Анализ кривых порогов слышимости показывает, что для построения аппроксимирующего уравнения необходимо использовать полином не ниже второго порядка. Варианты аппроксимации средних порогов слышимости уравнением при m 2: Из сравнения оценок показателя достоверности аппроксимации R 2 следует, что для формализованного описания кривых порогов слуха рассматриваемого пациента достаточно уравнения второго порядка m = 2: 52
Расчет коэффициентов динамических моделей На входе и выходе системы имеются зависимости X и Y от времени t Дискретизация: Таблица отсчетов, Δt = 0, 1 i 0 1 2 3 … i … n t 0 0, 1 0, 2 0, 3 … Δt · i … Δt · n xi 3 3, 2 3, 1 2, 6 … xi … xn Динамический ряд: совокупность значений переменной в таблице, упорядоченных во времени. Чем меньше расстояние между отсчетами, чем больше частота 53 дискретизации, тем меньше потери информации.
Расчет коэффициентов динамических моделей Любая динамическая система характеризуется коэффициентами производных (первой, второй и т. д. ) в записи модели. Чем большая степень старшей производной присутствует в записи модели, тем больший порядок динамической системы, тем глубже ее память, и тем больше коэффициентов (параметров) надо определить, чтобы идентифицировать систему. Оценить порядок динамической системы: совпадает со степенью наибольшей из производных Y по отношению к t. Система (звено) первого порядка При нулевых начальных условиях, если входной сигнал отсутствует, выходной сигнал равен нулю, система находится в покое. Если подать на вход единичный (пробный) сигнал и удерживать его на входе достаточно долго, то система на выходе начнет отклоняться от нулевого состояния. Система на выходе должна дойти до значения kx, k – коэффициент усиления входного сигнала. Насколько инерционно реагирует система, зависит от параметра T. Система достигнет значения kx на выходе и будет держать этот сигнал, пока на входе держится единичный сигнал. Переход от нуля до kx происходит во времени, – динамический процесс, т. е. в сигнале присутствует изменение, которое описывается производной, и выход оказывается меньше входа на некоторую величину: y = kx - f(dy/dt) 54
Расчет коэффициентов динамических моделей Звено первого порядка обладает двумя параметрами: инерционностью T и коэффициентом усиления k = Y(t = ∞)/X. T характеризует инерционность системы (память). При малой величине T система слабо зависит от предыстории и вход мгновенно заставляет измениться выход. При большом значении T система медленно реагирует на входной сигнал, а при очень большом значении T система выдает неизменный выходной сигнал, практически не реагируя на входные воздействия. Коэффициент k характеризует способность системы к усилению (при k < 1 – к ослаблению) уровня входного сигнала. Чтобы определить коэффициент k на графике, достаточно дождаться успокоения сигнала на выходе системы и вычислить отношение уровня выходного сигнала к уровню входного. Математически это означает, что все слагаемые, содержащие производные, равны нулю (система успокоилась, движения нет), а оставшееся слагаемое Y = k · X определяет значение k. 55
Расчет коэффициентов динамических моделей Звено первого порядка Передаточная функция как модель динамической системы, отношение выхода ко входу: W = Y/X Передаточная функция звена первого порядка: W = k / (Tp + 1), p – алгебраизованный оператор дифференцирования, символ дифференцирования, тождественно равный d/dt Y / X = k / (Tp + 1), откуда T · d. Y/dt + Y = k · X или T · ΔY/Δt + Y = k · X В разностном виде: T · (Yi + 1 – Yi) + Yi · Δt = k · Xi · Δt или, выразив настоящее через прошедшее: Yi + 1 = A · Xi + B · Yi, где A = k · Δt/T и B = 1 – Δt/T — весовые коэффициенты. A указывает на вес компоненты X, определяющей влияние внешнего мира на систему, B указывает на вес компоненты Y, определяющей память системы, влияние истории на ее поведение. Если B = 0, то Yi + 1 = А · Xi, безынерционная система Y = k · X, мгновенно реагирующая на входной сигнал и увеличивающая его в k раз. 56
Расчет коэффициентов динамических моделей Звено первого порядка Yi + 1 = A · Xi + B · Yi, A = k · Δt/T и B = 1 – Δt/T Если коэффициент B = 0, 5, т. е. 1 – Δt/T = 0, 5 или Δt/T = 0, 5, то A = k · Δt/T = k · 0, 5 и, следовательно, Yi + 1 = 0, 5 · k · Xi + 0, 5 · Yi. Экспонента, изображенная на графике, при большом n (в пределе n = ∞) стремится к значению входного (единичного) сигнала X, умноженного на коэффициент усиления k, что подтверждается расчетом: Yn + 1 = 0, 5 · k · Xn + 0, 5 · Yn = 0, 5 · k · Xn + 0, 5 · (0, 5 · k · Xn – 1 + 0, 5 · Yn – 1) = = … = (0, 51 + 0, 52 + … + 0, 5 n + 1) · k · X 0 + 0, 5 n + 1 · Y 0 = 1 · k · X 0. Выражение (0, 51 + 0, 52 + … + 0, 5 n + 1) является геометрической прогрессией, сумма которой при n = ∞ равна 1, стоящее при Y 0 выражение 0, 5 n + 1 обращается в 0 при n = ∞. 57
Расчет коэффициентов динамических моделей Звено второго порядка (колебательное звено) Чем больше производных учитывается в записи модели, тем больше порядок звена, тем больше коэффициентов при производных следует определить. Если на вход звена подать единичную функцию, при нулевых начальных условиях системы, то реакция на выходе будет называться переходной функцией (переходной характеристикой), h(t). Сигнал 1[t] – это, в некотором смысле, эталонный испытательный сигнал. Существуют и другие эталонные испытательные сигналы. Напр. , бесконечный импульс нулевой длины (дельта-функция Дирака), гармонический сигнал, периодические прямоугольные импульсы. Преобразуем уравнение: a 0 · p 2 · Y(p) + a 1 · p · Y(p) + a 2 · Y(p) = b · U(p) или: (a 0 · p 2 + a 1 · p + a 2) · Y(p) = b · U(p). Передаточная функция звена: T — постоянная времени (в секундах); ξ — коэффициент затухания (безразмерная величина); 58 k — передаточный коэффициент.
Расчет коэффициентов динамических моделей Звено второго порядка (колебательное звено) ξ ≥ 1 — апериодическое звено второго порядка 0 < ξ < 1 — колебательное звено второго порядка ξ = 0 — консервативное звено второго порядка 59
Расчет коэффициентов динамических моделей Динамические регрессионные модели, заданные в виде передаточной функции Примерный вид динамических сигналов на входе и на выходе; ограничимся временем рассмотрения сигналов T После дискретизации; отдельные отсчеты отстоят друг от друга на расстояние Δt, всего n отсчетов, т. е. T = n · Δt 60
Расчет коэффициентов динамических моделей Динамические регрессионные модели, заданные в виде передаточной функции Допустим, что передаточная функция W = Y/X (гипотеза должна быть подтверждена или опровергнута): или Дважды проинтегрируем это выражение и получим для некоторого произвольного момента времени t: Коэффициенты A 1, A 2, A 3, A 4 требуется определить. Выразим уравнение в разностном виде через суммы: где n – число экспериментальных точек. 61
Расчет коэффициентов динамических моделей Динамические регрессионные модели, заданные в виде передаточной функции Ошибка в некоторой m-ой точке: Ошибка показывает, насколько отходит теоретическое значение Ym от экспериментального значения. Суммарная ошибка минимизировать: (вносимая Величина ошибки зависит от значений параметров A 1, A 2, A 3, A 4, является функцией от четырех переменных. Чтобы найти минимум функции F, надо взять частные производные F по каждому из параметров и приравнять каждую производную к нулю. Система линейных уравнений: всеми точками), которую надо Из решения системы уравнений вычисляем неизвестные коэффициенты и дополняем ими модель: Необходимо: 1. сравнить получаемое из этой модели решение (Y теоретическое) с Y, заданным экспериментально, и вычислить ошибку F; 2. проверить, допустимо ли значение вычисленной ошибки, или гипотезу о виде модели 62 требуется сменить на более точную.
Моделирование биологических процессов Взаимодействие видов в биоценозах и взаимодействие химических веществ в растворах Пример: Модель химической реакции Лотки А. Д. (1926) Имеется химическая реакция, протекающая по общей схеме В некотором объеме находится в избытке вещество A. Молекулы A с некоторой постоянной скоростью k 0 превращаются в молекулы вещества X (реакция нулевого порядка). Вещество X может превращаться в вещество Y. Скорость этой реакции тем больше, чем больше концентрация вещества Y (зависит от продукта превращения, автокаталитический процесс). Молекулы Y необратимо распадаются, образуется вещество B (реакция первого порядка). Порядок реакции по данному веществу — показатель степени при концентрации этого вещества в кинетическом уравнении реакции. x, y, b – концентрации химических компонентов; k 0=ko’A, k 1, k 2 – константы скоростей реакций. Первые 2 уравнения этой системы не зависят от b, поэтому их можно рассматривать отдельно: 63
Моделирование биологических процессов Взаимодействие видов в биоценозах и взаимодействие химических веществ в растворах Пример: Экологическая модель Вольтерра В. (1931), модель хищник-жертва В некотором замкнутом районе живут жертвы и хищники, например зайцы и волки. Зайцы питаются растительной пищей, всегда имеющейся в достаточном количестве. Волки могут питаться лишь зайцами. Число зайцев – x, число волков – y. Т. к. количество пищи для зайцев не ограничено, можно предположить, что они размножаются со скоростью, пропорциональной их числу (аналогично автокаталитической реакции первого порядка): Пусть убыль численности зайцев пропорциональна произведению xy. Количество волков также нарастает тем быстрее, чем чаще их встречи с зайцами (тоже пропорционально xy, аналогично произведению концентраций). Кроме того, имеет место естественная смертность волков, причем скорость убывания особей пропорциональна их количеству. Тогда для изменений численности зайцев и волков получим: Кривые численности зайца и рыси в Канаде (К. Вилли, В. Детье, 1974) 64
Блок-схема автоматизированного расчета параметров модели Исходные данные расчета: - структура модели - начальные величины параметров - целевая функция - алгоритм минимизации целевой функции - воздействующие факторы 65
Оптимизация структуры модели Структура модели определяется множеством параметров (коэффициентов, переменных и отношений между ними). На заключительном этапе построения модели необходимо выбрать такую структуру, которая удовлетворяла бы заданным требованиям к модели, т. е. соответствовала бы некоторым критериям качества. При этом важно выбрать оптимальную по сложности модель. Существует 2 основных подхода к решению поставленной задачи: 1. Последовательное усложнение модели путем включения в нее новых переменных и перестройки взаимоотношений между ними. В случае, если новая переменная не играет существенной роли в повышении качества модели, то она заменяется новой переменной. Процесс заканчивается при достижении заданного качества модели. 2. Движение от сложного к простому, т. е. последовательное исключение малозначимых параметров, которые образуют «шум» в оценке расчетных значений зависимой переменной. Общим в этих подходах является то, что процедура оптимизации структуры модели представляет собой динамический процесс, в течение которого модель либо усложняется, либо упрощается. Оценка динамики критерия качества позволяет выбрать оптимальную структуру. Однако если критериев качества несколько, то один из критериев рассматривается как целевая функция, а для остальных задаются ограничения (стандартная задача оптимизации). 66
Оптимизация структуры модели Статистическая оценка значимости модели Необходимым условием для статистического оценивания модели является возможность сопоставления результатов эксперимента (с реальным объектом моделирования) и расчетных данных, полученных с помощью модели. Корреляция между расчетными и экспериментальными данными Оценка R 2 (коэффициент детерминации) Показывает, какая доля общей дисперсии зависимой переменной Y определяется расчетными (calc) значениями зависимой переменной, т. е. моделью Y = F(X): – математическое ожидание экспериментальных (exp) значений зависимой переменной; R – оценка коэффициента корреляции между расчетными и экспериментальными данными. Чем ближе R 2 приближается к единице, тем точнее построенная модель аппроксимирует экспериментальные данные. Оценка R 2 необходима для определения аппроксимируемости моделью экспериментальных 67 данных. Для медико-биологических исследований обычно необходимо R 2 ≥ 0, 5.
Статистическая оценка значимости модели Корреляция между расчетными и экспериментальными данными Линейный коэффициент корреляции указывает, есть ли между двумя рядами X и Y линейная зависимость и какой силы. mx, my, mxy — математическое ожидание x, y, xy: Дисперсия σx 2 и σy 2 показывает, насколько разбросаны точки относительно средней величины: Линейный коэффициент корреляции может иметь знак плюс или минус. Положительная его величина свидетельствует о прямой связи между X и Y. Чем ближе KR к +1, тем связь более тесная. Отрицательная величина его свидетельствует об обратной связи. Близость KR к нулю свидетельствует о слабой связи между X и Y. 68
Статистическая оценка значимости модели Информативность модели Оценивается для определения прогностической ценности модели. Применяется критерий Фишера: – расчетная дисперсия зависимой переменной; q – количество коэффициентов модели. – остаточная дисперсия, определяющая ошибку воспроизведения моделью реальных (экспериментальных) данных; N – количество опытов. Модель информативна, если из таблицы Фишера. Для входа в таблицу необходимо вычислить число степеней свободы: f 1 = q для расчетной дисперсии зависимой переменной, f 2 = N-q для остаточной дисперсии. Оценка информативности модели играет важную роль в тех случаях, когда нет возможности поставить дублирующие опыты при неизменных условиях (одинаковых величинах независимых переменных). Такая ситуация может возникнуть при моделировании зависимости какого-либо интегрального показателя состояния от комбинированного влияния ряда показателей 69 функционирования различных подсистем организма.
Статистическая оценка значимости модели Информативность модели Оценка информативности статистически определяет достоверность отличия модельных значений от «шума» , создаваемого неучтенными факторами. Оценки и информативности нередко вступают в противоречие: чем больше усложняется модель R 2 1, тем лучше она аппроксимирует экспериментальные данные, но тем меньше остается степеней свободы N-q и, следовательно, тем меньше вероятность того, что она информативна. Поэтому в процессе выбора структуры модели необходимо использовать методы оптимизации, обеспечивающие R 2 = max при Finf ≥ Finf допустимое. 70
Статистическая оценка значимости модели Адекватность модели В случае, когда имеется возможность оценить дисперсию воспроизводимости результатов эксперимента при повторном проведении опытов для одних и тех же условий, осуществляется проверка адекватности модели. Разброс экспериментальных точек относительно расчетных не должен превышать разброс в точках эксперимента относительно средних арифметических. Оценка разброса экспериментальных точек относительно расчетных – дисперсия адекватности: n – количество повторных опытов, N – количество опытов, – среднее арифметическое экспериментальных значений переменной в i-ом опыте, – расчетное значение переменной в i-ом опыте, q – количество коэффициентов модели. 71
Статистическая оценка значимости модели Адекватность модели Оценка разброса в точках эксперимента относительно средних арифметических (оценка точности эксперимента) – дисперсия воспроизводимости: – среднее арифметическое экспериментальных значений зависимой переменной в i-ом опыте, – значение переменной в i-ом опыте при j-ом дублировании измерения. Проверка адекватности модели производится на основе критерия Фишера: Гипотеза об адекватности модели принимается, если расчетное значение критерия Фишера не превышает его табличного значения при выбранном уровне значимости: Степени свободы: f 1 = N-q, f 2 = N(n-1) 72
Статистическая оценка значимости модели Адекватность модели Если опыты не дублируются (n = 1), то для оценки дисперсии воспроизводимости необходимо поставить один опыт с n повторными измерениями зависимой переменной: – среднее арифметическое результатов повторных измерений. Оценка дисперсии адекватности: Наиболее совершенной оценкой пригодности модели является проверка ее адекватности. Однако для этого необходимо проверить точность эксперимента, что далеко не всегда возможно, особенно при моделировании внутренних взаимоотношений в биологической системе. 73
Оценка значимости коэффициентов модели Для расчета коэффициентов экспериментально-статистических моделей с помощью МНК матричным способом* необходимо построить матрицу (XTX)-1. * B = (XTX)-1 XTY – коэффициенты модели Эта матрица составляет основу матрицы дисперсий-ковариаций (M-1), содержащей всю необходимую информацию для статистической оценки коэффициентов модели: Srep – дисперсия воспроизводимости; X – матрица значений независимых переменных; n – количество повторных опытов. Прямая матрица M называется информационной матрицей Фишера. Диагональные элементы матрицы – дисперсии коэффициентов, а недиагональные – ковариации между коэффициентами. Наличие ковариации свидетельствует о том, что оценки коэффициентов взаимосвязаны, что приводит к погрешности при их анализе и сопоставлении. Для однофакторной линейной полиномиальной модели y=b 0+b 1 x в общем виде матрица дисперсий-ковариаций: S 2(b 1) – дисперсия первого коэффициента модели, cov(b 0, b 1) – ковариация между коэффициентами b 0 и b 1. 74
Оценка значимости коэффициентов модели Важным условием построения матрицы дисперсий-ковариаций является невырожденность прямой матрицы (XTX), т. е. отсутствие нулей в знаменателе при ее обращении. Чтобы исключить это явление, необходимо отсутствие взаимной корреляции между независимыми переменными (X 1, X 2, …, Xm). Если при построении матрицы наблюдений X нет повторных (дублирующих) опытов, то вместо дисперсии воспроизводимости можно использовать дисперсию, характеризующую разброс экспериментальных значений зависимой переменной yexp относительно ее расчетных значений ycalc. Для любого коэффициента bi доверительный интервал: где t – табличное значение критерия Стьюдента (число степеней свободы берется то же, что и для дисперсии воспроизводимости), S(bi) – среднеквадратичная ошибка i-го коэффициента модели. Коэффициент bi значим, если его абсолютная величина больше доверительного интервала, т. е. Задача: Построить матрицу дисперсий-ковариаций и проанализировать значимость коэффициентов модели B при известных значениях этих коэффициентов, входных переменных X, дисперсии воспроизводимости и количестве повторных опытов. 75
Расчет ошибки прогнозирования Если Y = F(X 1, X 2) (линейная зависимость) то, согласно теории ошибок *, дисперсия расчетных значений зависимой переменной (ошибка предсказания, прогнозирования): * Прямая задача: перенос ошибок исходных данных (X) в результаты (Y) при произвольном преобразовании. В качестве переменных, влияющих на ошибку прогнозирования, можно рассматривать коэффициенты модели y = f(b 0, b 1, …, bi), тогда для простой модели y = b 0 + b 1 x: 76
Расчет ошибки прогнозирования Т. о. , оценка дисперсии предсказания Sy 2 зависит от трех основных факторов: - точности эксперимента; - координат точки предсказания (величин аргументов функции); - матрицы дисперсий-ковариаций, структура которой зависит от особенностей планирования эксперимента. Определение доверительных интервалов для предсказанного значения функции: где t – табличное значение критерия Стьюдента, определяемый из соответствующей таблицы; число степеней свободы то же, что и для дисперсии воспроизводимости. Задача: рассчитать ошибку прогнозирования и оценить доверительные интервалы выхода модели Y для всего диапазона входных переменных X при известных значениях входных переменных и матрице дисперсийковариаций. 77
Математические модели на основе активных экспериментов При исследовании многофакторных систем возможно 2 подхода: 1. Изучение влияния на выходную переменную каждого фактора в отдельности при фиксированных значениях остальных, всего (n!) комбинаций факторов. 2. Реализация специальных планов экспериментов, позволяющих выбрать оптимальное количество точек факторного пространства. В теории планирования экспериментов для описания зависимости выходной переменной y от факторов (x 1, x 2, …) используют понятие функции отклика. Геометрический образ, соответствующий этой функции, называется поверхностью отклика. Совокупность факторов образует так называемое факторное пространство. Поверхность отклика может быть представлена на факторной плоскости линиями постоянных значений функции отклика. Ф(x 1, x 2) x 1 x 2 78
Математические модели на основе активных экспериментов Планирование экспериментов связывают, прежде всего, с задачей минимизации числа необходимых экспериментов по оценке влияния на Y. При построении плана стремятся получить оптимальные свойства: - максимальную точность (минимальную дисперсию) в исследуемой области при заданном числе опытов, - независимость (некоррелируемость) оценок влияния разных факторов. Условия постановки опытов описываются в виде специальной таблицы, которая называется планом эксперимента. Для повышения эффективности экспериментальных исследований применяются специальные методы планирования, которые позволяют одновременно изменять все факторы и получать количественные оценки линейных эффектов и эффектов взаимодействия для уравнения регрессии. Существуют различные схемы планирования экспериментов, которые выбираются в зависимости от размерности модели, известного или предполагаемого распределения вероятностей ошибок измерения, вида поверхности отклика и т. д. 79
Полный факторный эксперимент (ПФЭ) предусматривает исследования влияния на выходную переменную изменения всех возможных комбинаций факторов на всех выбранных для исследования уровнях в некоторой локальной области. Необходимое число экспериментов для ПФЭ: где n – число факторов, L – число уровней варьирования. Уровнем фактора называют значение фактора, которое должно фиксироваться при проведении эксперимента. Наиболее простым является ПФЭ, в котором каждый фактор может принимать только два значения, L = 2. Обычно эти значения выбирают на границах диапазона изменения фактора. Т. о. , в ходе эксперимента каждый фактор принимает максимальное или минимальное значения. Полный перечень экспериментов с указанием значений факторов в каждом из них оформляется в виде специальной таблицы – матрицы планирования. Пример: t, o. C – температура P, атм. – давление X 1, X 2 – кодированные значения температуры и давления Y – значения зависимой переменной (скорости размножения бактерий) y 01 … m – результаты 80 параллельных экспериментов в центре плана
Полный факторный эксперимент В матрице планирования эксперимента каждая строка отвечает условиям проведения одного опыта, каждый столбец определяет значения одного фактора в разных опытах. План эксперимента можно рассматривать, как геометрическую интерпретацию расположения опытных точек в факторном пространстве. В полном факторном эксперименте рассматриваются все возможные комбинации факторов на всех выбранных для исследования уровнях. Обычно центр плана совмещают с началом координат факторного пространства. Эта операция эквивалентна преобразованию столбцов матрицы планирования в соответствии с выражениями: где Δzj – диапазон варьирования j-го фактора, zj 0 – значение j-го фактора в центре плана, zij – (натуральное) значение j-го фактора в i-ом опыте. Применив формулы ко 2 ому и 3 ему столбцам таблицы, получим кодированные значения факторов для ПФЭ. Координаты центра плана в кодированной матрице планирования с основанием 2 будут равны нулю. Обычно (при реализации линейной модели) матрицу ПФЭ дополняют кодированными значениями X 0, соответствующими коэффициентам свободного члена модели, 81 и равными +1.
Полный факторный эксперимент Матрица планирования с натуральными значениями факторов (столбцы 2, 3) используется для проведения экспериментов. Строки матрицы независимы и не регламентируют порядок проведения опытов. Очевидно, в таблице удобнее поменять местами строку 2 и строку 3 – сгруппировать опыты при одинаковой температуре, что позволит сократить общее время эксперимента. Матрица планирования может реализовываться m раз, тогда получится выборка параллельных значений y. С целью экономии общих затрат на экспериментальные исследования для постановки параллельных экспериментов и проверки воспроизводимости y можно использовать дополнительную точку плана, чаще всего для этого используют центр плана. Свойства матрицы планирования (в кодированном виде) Сумма элементов любого столбца, кроме нулевого и последнего (Y), равна нулю. Сумма квадратов элементов любого (кроме последнего) столбца всегда равна количеству экспериментов. Сумма произведений соответственных элементов двух любых столбцов (кроме последнего) равна нулю. => Свойство ортогональности: – единичная матрица, Ортогональность двух столбцов матрицы планирования означает полное отсутствие корреляции соответствующих факторов; следовательно, их оценки тоже не коррелированны. На основании экспериментальных данных и кодированных значений факторов оценки коэффициентов для линейной зависимости рассчитываются по формуле: 82
Полный факторный эксперимент Пример: Скорость размножения бактерий (штук в час) зависит от температуры и влажности воздуха. Требуется с помощью ПФЭ найти математическое описание этой зависимости в окрестности точки факторного пространства с координатами 50 о. С и 25 % с целью нейтрализации влияния этих факторов. Интервалы варьирования факторов: 5 о. С и 1 %. Чтобы убедиться в воспроизводимости процесса и уменьшить погрешность модели, каждый опытов выполняется дважды. Для каждой серии опытов вычисляется среднее значение и оценка дисперсии, которая определяет разброс значений функции отклика от среднего. Матрица планирования и результаты эксперимента: Номер опыта 1 2 3 4 1. 2. 3. 4. X 1 -1 +1 X 2 -1 -1 +1 +1 T, o. C 45 55 h, % 24 24 26 26 Yi 1 Yi 2 35, 0 39, 3 31, 8 36, 0 38, 1 33, 4 35, 6 Yiср 35, 5 38, 7 32, 6 36, 2 Дисп. 0, 5 0, 72 1, 28 0, 72 Определить коэффициенты регрессии (для кодированных значений факторов). Определить значимость коэффициентов регрессии. Проверить модель на адекватность. Получить модель в натуральном виде. 83
Дробный факторный эксперимент С увеличением числа факторов n резко возрастает число опытов в ПФЭ N. Если модель линейна, то для уменьшения объема экспериментальной работы можно использовать дробный факторный эксперимент (ДФЭ). План ДФЭ создаётся на основе определённой части плана ПФЭ, поэтому его иногда называют дробной репликой ПФЭ. Из множества n факторов ПФЭ отбирается (n–p) основных, и число необходимых опытов для постановки дробного факторного эксперимента: где p – показатель дробности, число линейных эффектов, приравненных к эффектам взаимодействия при построении плана ДФЭ (см. пример далее). Расчёт оценок коэффициентов регрессионной модели и проверка статистических гипотез по результатам ДФЭ осуществляются с использованием тех же соотношений, что и при постановке ПФЭ. Однако найденные значения b’i являются смешанными оценками генеральных коэффициентов регрессии βi. Особенности ДФЭ: - с уменьшением числа опытов N появляется корреляция между некоторыми столбцами матрицы планирования; - => становится невозможно раздельно оценивать эффекты факторов и эффекты взаимодействия; - => оценки b’i получаются смешанными. 84
Дробный факторный эксперимент Пример: Пусть надо построить математическую модель в виде линейного уравнения регрессии: Тогда для ПФЭ необходимо поставить 8 опытов N = 23, в ДФЭ число опытов можно уменьшить до четырех N = 23 -1 (полуреплика ПФЭ), т. е. 1 линейный эффект заменяется на эффект взаимодействия. Формирование матрицы ДФЭ: 2 ой и 3 ий столбцы матрицы планирования заполняются так же, как в ПФЭ, а 4 ый столбец формируется с помощью специального алгебраического соотношения – генерирующего соотношения (генератора плана). Генерирующим называется соотношение, показывающее, каким образом линейный эффект заменяется на эффект взаимодействия; оно задается как произведение основных факторов, определяющее значение элементов дополнительного столбца матрицы планирования. Для рассматриваемого примера возможно использование одного из генерирующих соотношений: С помощью столбцов 2– 4 и результатов эксперимента (выборки значений Y) можно найти оценки коэффициентов регрессии b’i, характеризующие линейные эффекты. Расчет осуществляется по тем же формулам, что и для ПФЭ. 85
Дробный факторный эксперимент Любая найденная на основе плана ДФЭ оценка коэффициента регрессии является смешанной. Правило смешивания, определяющее коррелированные основные эффекты и эффекты взаимодействия, описывают с помощью определяющего контраста реплики. Определяющий контраст полуреплики получается путем умножения генерирующего соотношения на его же левую часть, а т. к. для любой кодированной переменной xi 2=1, то левая часть формулы определяющего контраста всегда равна единице и обозначается I. В частности, для плана ДФЭ типа 23– 1 и генератора x 3 = x 1 x 2 имеет место определяющий контраст I = x 1 x 2 x 3 (генератор умножается на переменную x 3, следовательно, x 3 = I = x 1 x 2 x 3). Определяющий контраст позволяет найти, с какими параметрами смешана оценка коэффициента данного фактора. Для этого следует умножить обе части определяющего контраста на этот фактор. В рассматриваемой математической модели 3 фактора, поэтому получим 3 произведения: Оценки коэффициентов линейной модели будут смешанными: В рассматриваемом примере матрицу планирования для ДФЭ можно построить также с помощью генерирующего соотношения: Тогда смешанные оценки определяются соотношениями: 86
Дробный факторный эксперимент Реализовав обе полуреплики, путем совместной обработки результатов экспериментов можно получить раздельные оценки для линейных эффектов и эффектов взаимодействия (такой вариант плана соответствует ПФЭ): Разрешающая способность полуреплик (возможность раздельного определения коэффициентов модели) зависит от генерирующих соотношений. Так, если для плана 24– 1 выбрать генерирующее соотношение x 4 = x 1 x 2, то получим реплику с контрастом I = x 1 x 2 x 4 и разрешающей способностью x 1 = x 2 x 4 и т. д. Здесь линейные эффекты определяются совместно с парными взаимодействиями. Очевидно, что в первую очередь следует пренебречь взаимодействием более высоких порядков из-за их более низкой вероятности существования по сравнению с парными. Реплики можно строить высокой степени дробности, сокращая тем самым количество экспериментов. Взаимодействие факторов, выбранных в качестве генераторов плана, может быть значимым или незначимым. Для построения дробных реплик следует выбирать незначимые взаимодействия, которые определяются по физическим соображениям на основе априорных сведений. 87
Дробный факторный эксперимент Пример: Скорость размножения бактерий (штук в час) зависит от температуры, влажности воздуха и атмосферного давления. Требуется с помощью ДФЭ найти математическое описание этой зависимости в окрестности точки факторного пространства с координатами 50 о. С, 25 % и 760 мм. рт. с целью нейтрализации влияния этих факторов. Интервалы варьирования факторов: 5 о. С, 1 % и 10 мм. рт. ст. Чтобы убедиться в воспроизводимости процесса и уменьшить погрешность модели, каждый опытов выполняется трижды. Для каждой серии опытов вычисляется среднее значение и оценка дисперсии, которая определяет разброс значений функции отклика от среднего. Матрица планирования (полуреплика, генерирующее соотношение X 3= –X 1 X 2) и результаты эксперимента: Номер опыта X 1 X 2 X 3 T, o. C h, % p, мм. рт. ст. Yi 1 Yi 2 Yi 3 Yiср Дисп. 1 2 3 4 -1 +1 -1 -1 +1 +1 -1 45 55 24 24 26 26 750 770 750 105, 0 99, 3 115, 8 101, 8 106, 0 100, 1 115, 4 109, 8 99, 9 115, 1 110, 6 103, 6 99, 8 15, 4 108, 6 11, 08 0, 173 0, 123 7, 213 1. 2. 3. 4. Определить коэффициенты регрессии (для кодированных значений факторов). Определить значимость коэффициентов регрессии. Проверить модель на адекватность. 88 Получить модель в натуральном виде.
Построение моделей на основе планов второго порядка Уравнения математических моделей, полученных путем полного или дробного факторных экспериментов используют для решения двух типов задач: прогнозирования и оптимизации. Задача прогнозирования заключается в определении оценки выходной переменной для заданных значений факторов. Очевидно, если найдены оценки коэффициентов уравнения модели, значение Y получается простой подстановкой значений факторов в уравнение. Задача оптимизации сводится к поиску таких значений факторов, при которых Y или некоторый функционал F(Y) принимают экстремальное значение. При решении задач оптимизации осуществляется направленное движение в факторном пространстве от некоторой исходной точки по направлению к точке экстремума. Пример задачи оптимизации: Исследования факторного пространства вблизи точки экстремума целевой функции показало, что уравнение регрессии, полученное на основании ПФЭ или ДФЭ, не позволяет адекватно описывать объект. Модель, адекватная в области вдали от экстремума целевой функции, не обеспечивает удовлетворительного по точности прогноза в области ее экстремума. В окрестностях экстремума целевой функции довольно сильно проявляются эффекты взаимодействия и квадратичные эффекты. Для описания Y в этой области используют полиномы второго порядка: 89
Построение моделей на основе планов второго порядка Для определения оценок коэффициентов регрессии в таком случае необходимы эксперименты, в которых каждый фактор варьировался бы не менее, чем на трех уровнях, L ≥ 3. Так как исследователь обращается к планам второго порядка обычно после того, как не удалось получить адекватную модель на основе ПФЭ (или ДФЭ), то естественно желание сохранить результаты, полученные в этих экспериментах, и использовать их в дальнейших расчетах. С учетом этих соображений разработаны композиционные планы 2 -го порядка. 90
Построение моделей на основе планов второго порядка Рассмотрим построение регрессионной модели на основе ортогонального центрального композиционного плана (ОЦКП). Ядром этого плана служат точки ПФЭ или ДФЭ. Центр ОЦКП совпадает с центром плана, принятого в качестве ядра. Новый план дополняется так называемыми звездными точками, которые располагаются на координатных осях на расстоянии звездного плеча α от центра плана. В состав плана добавляется также одна или несколько точек в центре плана. Расположение точек в ОЦКП при двух независимых факторах: Для двухфакторного эксперимента точки плана лежат в вершинах квадрата, вписанного в окружность – проекцию параболоида на плоскость факторов. Для задания линейной модели, которая представляет собой уравнение прямой или семейства прямых достаточно задаваться лишь двумя уровнями факторов (+1 и – 1), т. к. прямую можно однозначно провести через две известные точки. В модели 2 -го порядка поверхность отклика образуется семейством кривых 2 -го порядка, для построения которых необходимо знать не менее 3 -х точек. Следовательно, в данном случае переменные нужно варьировать не менее чем на 3 -х уровнях. - 91
Построение моделей на основе планов второго порядка Общее число точек плана второго порядка определяется по формуле где Nс – число точек ядра плана (точки плана ПФЭ), Nst – число звездных точек, N 0 – число опытов в центре плана ядра. Рассмотрим описание выхода модели Y в виде полинома второго порядка: Построим матрицу композиционного плана для двух факторов, в которой предусматривается единственный опыт в центре плана, n = 2, N 0 = 1. Используя в качестве ядра матрицу ПФЭ, получим неортогональную матрицу: Условия ортогональности нарушаются для столбцов: Произведя замену квадратичных переменных и выбрав величину плеча α, можно преобразовать матрицу плана к ортогональному виду. 92
Построение моделей на основе планов второго порядка Приведем матрицу к ортогональному виду. Для этого преобразуем квадратичные переменные по формуле Это эквивалентно введению новых переменных и ведет к изменению условий ортогональности: Из последнего условия получают уравнение для звездного плеча: Для n = 2, N 0 = 1, = 1 Для определения звездного плеча можно пользоваться готовыми таблицами, составленными для разного числа факторов. На число опытов в центре плана ограничений не накладывается. В силу ортогональности модифицированной матрицы независимых переменных все коэффициенты уравнения модели определяются независимо один от другого. Затем уравнение модели пересчитывается для старых переменных. 93
Построение моделей на основе планов второго порядка Для определения коэффициентов регрессии в квадратичной модели удобно использовать композиционный (составной) план, состоящий из: • полной или дробной реплики симметричного двухуровневого плана; • центральной точки, расположенной в начале координат факторного пространства; • звёздных точек. № опыта 1 2 3 4 5 6 7 8 9 X 1 +1 -1 0 0 +1, 58 -1, 58 0 X 2 +1 +1 -1 -1 +1, 58 -1, 58 0 0 0 Реплика ПФЭ Звёздные (. . ) Центральная (. ) 94
Моделирование биологической системы Пример Постановка задачи: Необходимо построить модель газообмена в легких человека. Формулировка задачи: «В каком количестве присоединяется гемоглобином кислород в крови в зависимости от внешних влияющих факторов? » Известно, что для эффективного газообмена в легких необходимо поддержание определенного трансмурального давления (разности давлений в сосуде и окружающей среде) в капиллярах кровеносной системы. ОИ как система: Выход модели – распределение трансмурального давления по сосудам кровеносной системы. Процесс дыхания не сводится к движению воздуха в бронхолегочной системе, а протекает при постоянном газообмене в легочной ткани, что обеспечивается интенсивным обмыванием легких кровью. При этом давление крови в капиллярах оказывает существенное влияние на эластичные свойства легочной ткани. => Для полного и точного описания дыхательных процессов в живом организме необходимо математически описать процессы, протекающие в кровеносной системе. Сосудистая система (в частности малого круга кровообращения) является не просто системой трубопроводов с фиксированными сечениями каналов, определяемых «жесткой» морфометрической структурой, а является автономной саморегулируемой гидродинамической системой стабилизации распределения давления вдоль всего тракта. 95
Модель газообмена в легких человека Построение модели: Идеализируя структуру объекта исследования и ближайшего функционального пространства, представим их в виде модели, в которой: • легкие имеют вид сосуда переменного объема; • внутри объема легких в упругой среде размещена морфометрическая структура кровеносных сосудов. Аппроксимационные зависимости длин и диаметров сосудов от порядка ветвлений: для артериальных сосудов для венозных сосудов • кровеносная система связана через предсердие с сердцем, выполняющим функцию насоса, прокачивающего кровь по морфометрической структуре. Малый круг кровообращения начинается легочным стволом, который является продолжением артериального конуса правого желудочка сердца, разделяется на правую и левую легочную артерии, внедряющиеся в ткань легких в области их корней. Ветвление легочной артерии происходит в соответствии с ветвлениями бронхиального дерева вплоть до респираторных бронхиол, на уровне которых формируются мельчайшие артерии и артериолы для каждого альвеолярного хода, заканчивающиеся сетью капилляров на поверхности альвеолярных мешочков и альвеол. Описание артериального отдела сосудистой системы малого круга предусматривает его деление на 17 генераций. В общем случае артериальное дерево легких описывается морфометрической моделью неправильного ветвления. При этом давление и скорость движущейся жидкости (крови) в системе сосудов является неоднозначной функцией порядка ветвления сосудов => модели правильного ветвления (более простая). 96
Модель газообмена в легких человека Приведенная схема замещения включает окружающее объект моделирования функциональное пространство. Граничные условия выделения системы: 1. Деформация легочных объемов связана с работой сил, создаваемых мышечной структурой легких. Эти силы оказывают воздействие на окружающую кровеносную систему. Если задать распределенное поле давлений, непосредственно воздействующее на отдельные сосуды как воздействие легкого учитывается как опосредованное (через заданное распределение давлений). 2. Основная функция сердца в рассматриваемой схеме замещения – насосная, обеспечивающая средний объемный расход крови за определенный интервал жизнедеятельности организма, определяемый потребностью в химических веществах: 3. Для обеспечения соответствующего расхода крови в легочной артерии создается давление, превышающее давление в венозном русле. Если в качестве переменных модели задать расход, артериальное давление и венозное давление как функции времени, то система малого круга кровообращения (как гемодинамическая система) может рассматриваться отдельно от окружающего ее функционального пространства. Подсистемы: артерии эластичного типа, служащие поддерживающими структурами; артерии мышечного типа, обладающие свойствами сократимости; артериолы и венулы, обладающие 97 упругими свойствами и деформирующиеся в процессе дыхания. Свойства сосудов: упругость (эластичность) и мышечная сократимость.
Модель газообмена в легких человека Предложенная модель не учитывает тот факт, что мы имеем дело с гидродинамической системой, в которой кровь является «тяжелой жидкостью» и, находясь в поле тяготения, будет изменять гидростатическое давление. Таким образом, на различных участках малого круга кровообращения в легких под действием поля тяготения будет возникать различное гидростатическое давление (относительно уровня сокращения сердца). Для учета гравитационного воздействия схема замещения должна учитывать расположение отдельных участков морфометрической структуры в легких относительно уровня сокращения сердечной мышцы. 98
Модель газообмена в легких человека При течении жидкости в системе сосудов с неправильным ветвлением давление и поток не являются однозначными функциями порядка ветвления сосудов => модель правильного ветвления: каждый сосуд разветвляется на k одинаковых сосудов, имеющих на единицу больший номер поколения. Все сосуды одного поколения имеют одинаковые размеры и ответвляются под одинаковыми углами от отцовских сосудов. Число k называется отношением ветвления. Модель правильного ветвления считается заданной, если задано отношение ветвления и зависимость длины и диаметра сосудов от номера поколения. В системе сосудов с правильным ветвлением номер поколения сосуда однозначно определяет его размеры и положение по отношению к началу ветвления дерева сосудов, а поток и давление являются однозначными функциями номера поколения сосуда. Модель системы легочных сосудов – структура с правильной k-хотомией, т. е. с числом ветвления k. На каждом j-ом уровне ветвления кровеносного русла имеется Nj одинаковых участков: где n – количество рассматриваемых уровней морфометрической структуры. На j-ом уровне ветвления во всех Nj участках параметры движения жидкости (скорость, давление, трение) одинаковы. 99
Модель газообмена в легких человека Полагая, что средняя скорость движения жидкости в одном j-ом участке равна скорости движения центра масс объема, заключенного в данном участке: (1) где Sj – площадь поперечного сечения одного j-ого участка, mj = Sj lj j – масса жидкости в объеме одного j-ого участка, j – плотность жидкости на j-ом участке, Vj – скорость жидкости на j-ом участке, pj и pj+1 – статические давления соответственно в начале и конце участка, zj и zj+1 – высота соответственно начала и конца участка, отсчитанная от некоторого фиксированного уровня, Fтрj – сила вязкого трения. 100
Модель газообмена в легких человека Сила вязкого трения должна учитывать как потери на трение жидкости о стенки сосудов, так и местные потери при ветвлении. Считая, что разветвление потока жидкости происходит плавно, без внезапного расширения или сужения, потерями по полному давлению от местных гидродинамических сопротивлений можно пренебречь, и тогда сила трения будет определяться только вязким трением жидкости о стенку: – коэффициент сопротивления, зависящий от числа Рейнольдса Re = VD/ , , – динамический и кинематический коэффициенты вязкости, = . 101
Модель газообмена в легких человека Зависимость вязкости j от номера участка ветвления j: где – кинематический коэффициент вязкости крови в крупных сосудах, пределяемый экспериментальным путем, – коэффициент снижения вязкости крови, , – аппроксимационные коэффициенты. Как показывают расчеты, движение крови в сосудах отвечает условию ламинарного течения для которого коэффициент и сила трения определяется как Тогда уравнение (1) после элементарных преобразований приводится к виду: Вводя статические давления на входе и выходе морфометрической структуры pвх и pвых соответственно, высоты расположения входа zвх и выхода zвых, проводя суммирование по количеству уровней и считая кровь несжимаемой жидкостью с постоянной плотностью, получим: 102
Модель газообмена в легких человека Для осуществления дальнейших преобразований уравнения движения рассмотрим закон сохранения массы рабочего тела на каждом j-ом уровне морфометрической системы. При этом изменение массы рабочего тела на каждом j-ом уровне определим как разность прихода массы из предшествующего и ухода из рассматриваемого уровня, выраженные через средние скорости (скорости движения центров масс) движения жидкости на данных уровнях, т. е. (2) где mcj = Nj j Sj lj – масса жидкости во всех сосудах j-ого уровня. Полагая, что все процессы, связанные с изменением геометрии каналов и плотности газа в них, являются внешними и могут рассматриваться по отношению к процессу движения как квазистационарные, принимаем Тогда из уравнения (2) следует, что величина Nj j. Sj. Vj одинакова на всех уровнях и равна массовому расходу G рабочего тела через морфометрическую систему, т. е. G = Nj j Sj. Vj. Разрешая последнее уравнение относительно комплекса ( j Sj Vj) получим: (3) 103
Модель газообмена в легких человека Движение крови в сосудах легких будем рассматривать как движение несжимаемой жидкости из объема в объем. Один из объемов - сердечная камера, второй объем – капиллярная сеть малого круга кровообращения. В выделенных объемах жидкость практически малоподвижна, т. е. статические давления в объемах практически равны полным за исключением участков вблизи входных (выходных) отверстий – клапанов сердца и капиляров. Статическое давление при истечении жидкости из объема определяется из уравнения Бернулли для потока несжимаемой среды в этом сечении: и условия, что полное давление во входном сечении равно полному давлению в объеме (p 0 вх = p 0 w). При втекании жидкости в объем статическое давление в выходном сечении канала равно статическому давлению в объеме и, следовательно, полному давлению в объеме (pвых = pw 2 = p 0 w 2). Вводя величины pвн 1 и pвн 2 – внутренние давления в объемах, уравновешивающие внешние давления на жидкость в этих объемах, в частности, давления, вызываемые сокращением мышц сердца и альвеолярное давление на капиллярные сосуды, получим: (4) Подставляя уравнение (4) в уравнение (3) и, учитывая постоянство плотности крови, переходя к объемному расходу Q = G/ , будем иметь: где введено обозначение p = pвн 1 – pвн 2 + g (zвх – zвых). (5) 104
Модель газообмена в легких человека Уравнение (5), записанное для артериальной и венозной ветвящейся сети кровеносных сосудов при заданных внешних давлениях на входе и выходе (в сердце и в легких) полностью описывает движение крови, в том числе и неустановившееся, по кровеносной системе в целом. При этом динамические параметры течения - скорости и статические давления на каждом j-ом участке не входят в эти уравнения в явном виде. Если необходимо определить эти параметры на каком - либо участке, например k-ом, то можно поступить следующим образом. Проводя суммирование не по всем участкам, а только по участкам с номерами 1. . . (k-1), можно получить а после преобразований, аналогичных проводимым при выводе уравнения (5) можно прийти к уравнению откуда (6) Это уравнение позволяет определить статические давления на любом участке ветвления кровеносной системы. При этом значение объемного расхода Q получается при интегрировании дифференциального уравнения (5), а величина d. Q/dt, входящая в уравнение (6), подставляется из уравнения (5) непосредственно. 105
Методы исследования моделей Модель рассматривается как средство или «инструмент» для изучения ОИ и решения поставленных задач. Методы исследования моделей должны быть направлены на получение необходимой информации о моделируемом объекте и на принятие решений по результатам анализа. В процессе исследования модели необходимо: - выбрать методы расчета и анализа выходных параметров; - наметить стратегию проведения экспериментов с моделью. Методы расчета и анализа выходных параметров моделей Для математических моделей: - аналитические; - качественные; - численные. Для имитационных моделей: проведение машинных экспериментов (актуальна задача оптимизации проведения экспериментальных исследований имитационных моделей). 106
Качественные методы исследования моделей позволяют оценить характеристики и поведение системы качественно, минуя аналитическое решение математического выражения, которое ее описывает; используются для исследования динамических систем, представленных в виде ДУ, и статических систем, моделируемых с помощью многомерных нелинейных полиномов. Динамические системы При моделировании биомедицинских систем, когда значения параметров и условий часто не могут быть точно заданы, определяют их существенные качественные характеристики: существуют ли в системе стационарные состояния равновесия; устойчивы ли они; как меняется характер устойчивости в зависимости от параметров. Метод фазовой плоскости Рассмотрим систему ОДУ, где P, Q – функции, непрерывно дифференцируемые в некоторой области D изменения переменных x 1 и x 2, t – время: (1) 107
Метод фазовой плоскости Системе присущи 3 типа фазовых траекторий на фазовой плоскости x 1 Ox 2: точка, замкнутая кривая и незамкнутая кривая. Фазовая плоскость – это двумерное пространство, координатами которого являются величины неизвестных функций. Точкам соответствуют положения равновесия – такие решения системы, при которых фазовые переменные x 1, x 2 не изменяются во времени t, т. е. Замкнутые кривые на фазовой плоскости изображают решения, описываемые периодическими функциями (периодические решения) с периодами T 1 и T 2 соответственно, удовлетворяющими условию m. T 1=n. T 2 (m, n – натуральные числа). Незамкнутые кривые на фазовой плоскости изображают решения системы ДУ (1), которые не являются периодическими. Основная задача при качественном исследовании системы – установление структуры разбиения фазовой плоскости на траектории. Фазовый портрет системы – картина, которую образуют фазовые траектории на плоскости x 1 Ox 2. 108
Метод фазовой плоскости Каждая из фазовых траекторий является годографом вектора (2) решения системы (1), если начало этого вектора зафиксировано в начале системы координат x 1 Ox 2, т. е. если это радиус-вектор. Тогда x(t) – фазовый вектор. Годограф – кривая, представляющая собой геометрическое место концов переменного вектора, значения которого в разные моменты времени отложены от общего начала. Рассмотрим: Однородная система линейных ОДУ с постоянными коэффициентами: (3) Если система имеет единственное положение равновесия в начале координат: (4) Этой точке фазовой плоскости отвечает тривиальное решение системы (3): которое можно рассматривать как невозмущенное движение. 109
Метод фазовой плоскости Характер поведения фазовых траекторий системы (3) на фазовой плоскости x 1 Ox 2 определяют корни λ 1 и λ 2 характеристического уравнения: (5) Эти корни могут быть либо действительными причем простыми или кратными либо комплексно сопряженными 1. Корни действительные простые В этом случае общее решение системы (3) в векторной форме имеет вид (6) где C 1 и С 2 – постоянные, a 1 и a 2 – собственные векторы матрицы системы ОДУ. Векторы a 1 и a 2 линейно независимы, поскольку отвечают различным собственным значениям λ 1 и λ 2 этой матрицы и, следовательно, образуют базис на фазовой плоскости. Тогда (7) являются координатными функциями вектор-функции x(t) в этом базисе. В общем случае направления векторов a 1 и a 2 не совпадают с направлениями координатных осей Ox 1 и Ox 2 и система координат ξ 1 O ξ 2, вообще говоря, не является прямоугольной. Однако для наглядности удобно использовать переменные ξ 1, ξ 2 как прямоугольные координаты другой плоскости и изображать поведение 110 системы на этой плоскости в прямоугольной системе координат ξ 1 O ξ 2.
Метод фазовой плоскости В системе координат ξ 1 O ξ 2 фазовую траекторию в параметрической форме задают соотношения (7), причем роль параметра играет время t. Если положение точки на фазовой траектории характеризовать значением t, то при помощи (7) можно проследить за ее движением по фазовой траектории, при этом движущуюся точку называют изображающей. Ясно, что координаты x 1, x 2 и ξ 1, ξ 2 одной и той же точки фазовой траектории связаны линейным преобразованием, которое несколько видоизменяет фазовый портрет при переходе от одной системы координат к другой, но сохраняет его основные особенности. Эти особенности существенно зависят от знаков корней λ 1 и λ 2. 1. 1. Корни отрицательны (7) => при т. е. изображающая точка по любой фазовой траектории стремится к началу координат, но начало координат формально принадлежит лишь одной фазовой траектории (при C 1=C 2=0), которая вырождается в точку и соответствует положению равновесия (4). Если то согласно (6), (7) в системе координат ξ 1 O ξ 2 фазовая траектория – полуось оси O ξ 1, а если наоборот, то оси O ξ 2. 111
Метод фазовой плоскости Любой фазовой траектории, задаваемой соотношениями (7), в первом квадранте плоскости ξ 1 O ξ 2 при конкретных значениях отвечают еще 3 фазовые траектории: симметричная относительно оси O ξ 1 при смене знака постоянной C 2, симметричная относительно оси O ξ 2 при смене знака постоянной C 1 и симметричная относительно начала координат при смене знаков C 1 и C 2 одновременно. Следовательно, чтобы построить фазовый портрет системы (3) в плоскости ξ 1 O ξ 2, достаточно построить ее фазовые траектории в первом квадранте этой плоскости. Исключив при время t, получим Т. о. , фазовые траектории в плоскости ξ 1 O ξ 2 имеют вид ветвей парабол. Если к ним присоединить начало координат, то в нем они будут иметь общую касательную – ось O ξ 1. Фазовый портрет в переменных ξ 1, ξ 2 и x 1, x 2 112
Метод фазовой плоскости При отрицательных значениях корней λ 1 и λ 2 линейная система ОДУ асимптотически устойчива, т. к. все решения этой системы, асимптотически устойчивы, что позволяет говорить об устойчивости положения равновесия системы. На фазовой плоскости оно отвечает началу координат, куда по каждой из фазовых траекторий стремится изображающая точка при В данном случае положение равновесия называют устойчивым узлом. 1. 2. Корни положительны В этом случае, рассуждая аналогично случаю отрицательных корней, придем к выводу, что фазовые траектории имеют ту же форму. Однако при движение изображающей точки по фазовой траектории происходит от начала координат. Поэтому фазовый портрет совпадает с рис. на слайде 112, но стрелки направлены в обратном направлении. Решение системы (3), а значит и положение равновесия в начале координат, при положительных корнях будет неустойчивым по Ляпунову. В этом случае говорят, что положение равновесия в начале координат, называемое неустойчивым узлом, не является устойчивым. 113
Метод фазовой плоскости 1. 3. Корни имеют разные знаки В случае получим т. е. фазовая траектория – полуось оси O ξ 2, а движение изображающей точки происходит от начала координат. Если то т. е. фазовая траектория – полуось оси O ξ 1, а движение изображающей точки происходит к началу координат. При уравнение фазовой траектории, полученное из (7) исключением времени, будет Т. е. фазовые траектории имеют вид гипербол, для которых оси O ξ 1, O ξ 2 являются асимптотами. Положение равновесия в начале координат, соответствующее невозмущенному движению – седло. Седлу соответствует невозмущенное движение, которое не будет асимптотически устойчивым решением системы и не является устойчивым положением равновесия. 114
Метод фазовой плоскости 1. 4. Один из корней равен нулю (7) => Это тождество задает в плоскости ξ 1 O ξ 2 семейство прямолинейных фазовых траекторий, параллельных оси O ξ 1. Направление движения по ним изображающей точки зависит от знака корня λ 1. Если то ее движение происходит в направлении к оси O ξ 2. Если то в противоположном. Из характеристического уравнения (5) следует, что один из корней равен нулю при условии т. е. когда матрица A системы ОДУ является вырожденной. (8) В этом случае однородная система линейных алгебраических уравнений имеет бесконечное множество решений, которым соответствует бесконечное множество положений равновесия. Действительно, всякая точка оси O ξ 1, т. е. прямой ξ 1=0, является в данном случае положением равновесия системы (8), и отвечающие им невозмущенные движения при являются устойчивыми по Ляпунову, но не будут асимптотически устойчивыми, а при будут неустойчивыми по Ляпунову. 115
Метод фазовой плоскости 2. Корни действительные кратные В этом случае могут возникнуть 2 различных варианта в зависимости от значений коэффициентов a 12 и a 21 в правой части (3). 2. 1. При система является совокупностью 2 ух независимых линейных однородных ОДУ первого порядка с решениями: (9) Тогда исключив из (9) время, получим уравнение пучка прямых в плоскости x 1 Ox 2, проходящих через начало координат, которое следует рассматривать как положение равновесия (4) такой системы. Каждая полупрямая, не содержащая начало координат, будет фазовой траекторией, а изображающая точка будет приближаться по ней к началу координат, если и удаляться, если Положение равновесия в начале координат называют дикритическим узлом. Такое положение равновесия является асимптотически устойчивым при и не будет устойчивым при Частный случай возможен, если коэффициенты первой части системы (3) равны нулю. Тогда любая точка фазовой плоскости является положением равновесия такой системы, а соответствующее ему невозмущенное движение устойчиво по Ляпунову, но не является асимптотически устойчивым. 116
Метод фазовой плоскости 2. 2. Пусть хотя бы один из коэффициентов a 12 и a 21 в правой части (3) не равен нулю. Исключим из системы (3) переменную x 2, для это продифференцируем первое уравнение системы: Подставляя в это соотношение вместо dx 2/dt правую часть второго уравнения системы (3), получим: Выразив a 12 x 2 через x 1 при помощи 1 го уравнения системы (3), получим линейное однородное ОДУ второго порядка с постоянными коэффициентами, (10) характеристическое уравнение которого, естественно совпадает с (5) и имеет в рассматриваемом случае кратные корни Общее решение ОДУ (10) будет (11) Подставив (11) в первое ОДУ системы (3), найдем (12) 117
Метод фазовой плоскости Используя (11) и (12), общее решение системы (3) в случае корня λ кратности 2 запишем в виде (13) где – координатные функции векторфункции x(t) в базисе, образованном линейно независимыми векторами При имеем положение равновесия системы в начале координат. Если то фазовые траектории являются полуосями оси Oξ 1. При имеем и в системе координат ξ 1 Oξ 2 получаем уравнение фазовой траектории в виде: В случае т. е. изображающая точка по каждой из фазовых траекторий неограниченно приближается к положению равновесия (4) в начале координат, но не достигает его. Фазовый портрет для переменных ξ 1, ξ 2 и x 1, x 2 (при произвольной ориентации линейно независимых векторов a 1 и a 2) 118
Метод фазовой плоскости В случае фазовый портрет аналогичен, но направление движения изображающей точки противоположно. Положение равновесия в обоих случаях называют вырожденным узлом. Оно асимптотически устойчиво при и не является устойчивым при В случае координатные функции принимают вид Прямолинейные фазовые траектории параллельны оси Oξ 1, все точки которой являются положениями равновесия. Направление движения изображающей точки по фазовой траектории зависит от знака C 2. 3. Корни комплексно сопряженные Условием возникновения комплексно сопряженных корней характеристического уравнения (5) является выполнение неравенства или Отсюда следует, что коэффициенты a 12 и a 21 должны быть отличны от нуля и иметь разные знаки. Тогда 119
Метод фазовой плоскости Действительное решение ОДУ (10), соответствующего системе (3), можно записать в виде (14) Подставив (14) в первое ОДУ системы (3), получим (15) Используя (14) и (15), общее решение системы (3) в случае комплексно сопряженных корней запишем в виде (13), но теперь координатными функциями вектор-функции x(t) в базисе линейно независимых векторов являются (16) Вместо (16) можно написать (17) где Рассмотрим возможные частные случаи. 120
Метод фазовой плоскости 3. 1. Если то и корни характеристического уравнения являются чисто мнимыми числами Исключая в этом случае из (17) время t, получаем уравнение фазовой траектории Т. о. , в плоскости ξ 1 Oξ 2 фазовые траектории являются семейством концентрических окружностей с центром в начале координат, которые в плоскости x 1 Ox 2 переходят в семейство эллипсов, т. е. фазовые траектории замкнуты, что свидетельствует о том, что решение системы периодическое. При выбранных координатных функциях ξ 1(t), ξ 2(t) и ориентации координатных осей в плоскости ξ 1 Oξ 2 направление движения изображающей точки происходит против часовой стрелки. В плоскости x 1 Ox 2 это направление зависит от взаимного расположения векторов a 1 и a 2. Начало координат соответствует значениям C 1=C 2=0 (ρ=0) и является положением равновесия системы. Его в этом случае называю центром. Оно устойчиво, 121 хотя и не асимптотически.
Метод фазовой плоскости 3. 2. Если комплексно сопряженные корни имеют отрицательную действительную часть тогда соотношения (17) будут в параметрической форме задавать фазовые траектории в виде логарифмических спиралей, по которым изображающая точка будет стремиться к началу координат при В данном случае начало координат является асимптотически устойчивым положением равновесия системы и его называют устойчивым фокусом. 3. 3. Если комплексно сопряженные корни имеют положительную действительную часть то фазовый портрет аналогичен, но направление движения изображающей точки по фазовой траектории противоположно, поскольку при Положение равновесия в этом случае не является устойчивым и его называют 122 неустойчивым фокусом.
Метод фазовой плоскости 123
Метод фазовой плоскости Построение фазового портрета методом изоклин Если в системе ДУ (1) второе уравнение разделить на первое, получим уравнение, которое в явном виде не содержит время: dx 2/dx 1=Q/P. Изоклины – линии на фазовой плоскости, которые пересекаются со всеми фазовыми траекториями од одним и тем же углом, т. е. dx 2/dx 1=const. Особый интерес представляют главные изоклины, для которых: dx 2/dx 1=0 и dx 1/dx 2=0 (или dx 2/dx 1=∞) или Q(x 1, x 2)=0 и P(x 1, x 2)=0. Первая из этих изоклин пересекается с фазовыми траекториями параллельно оси x 1, а вторая – параллельно оси x 2. Пересечение главных изоклин образует особую точку, которая соответствует состоянию равновесия системы. Направления пересечения главных изоклин определяются из неравенств: Q(x 1, x 2)>0 (движение вверх) и P(x 1, x 2)>0 (движение вправо). 124
Метод фазовой плоскости Модель хищник-жертва Исследуем особую точку вольтерровской модели хищник-жертва. Ее координаты легко найти, приравняв к нулю правые части системы ДУ, описывающей модель. Это дает стационарные ненулевые значения: Т. к. все параметры положительны, точка будет расположена в положительном квадранте фазовой плоскости. Линеаризация системы вблизи этой точки дает: отклонения от особой точки на фазовой плоскости. Характеристическое уравнение системы Корни чисто мнимые Фазовые траектории вблизи особой точки представляют собой концентрические эллипсы. В целом, особая точка не является устойчивой: Колебания системы изменяют свои характеристики при изменении внешнего воздействия. 125
Численные методы решения дифференциальных уравнений Моделирование БТС => проблемы, связанные с решением ДУ аналитическими методами (нелинейность уравнений, их высокий порядок, наличие переменных коэффициентов). Численные методы: элементарные арифметические операции; частные решения для заданных начальных условий в табличном виде (дискретные последовательности, а не непрерывные кривые). Для решения ДУ численными методами их необходимо преобразовать в алгебраические уравнения в конечных разностях => конечно-разностные методы. Метод Эйлера Пусть нам известна входная динамическая последовательность X (входной сигнал) и модель (способ преобразования входного сигнала в выходной сигнал). Рассматривается задача определения выходного сигнала y(t). Модель динамической системы может быть представлена дифференциальным уравнением. Основное уравнение динамики: y' = f(x(t), y(t), t). Известны начальные условия в нулевой момент времени t 0: y(t 0), x(t 0). По определению производной: 126
Метод Эйлера Иллюстрация расчета будущего состояния системы методом Эйлера на одном шаге Нам известно положение системы в точке « 1» , требуется определить положение системы в точке « 2» . Точки отделены друг от друга расстоянием Δt. То есть расчет поведения системы производится по шагам. Из точки « 1» мы скачком (дискретно) переходим в точку « 2» , расстояние между точками по оси t называется шагом расчета Δt. Тогда: или: Формула Эйлера 127
Метод Эйлера Геометрическая иллюстрация метода Эйлера Пусть A — точка, в которой состояние системы известно. Это «настоящее» состояние системы. В точке A к траектории движения системы проведем касательную. Касательная — это производная функции f(x(t), y(t), t) по переменной t. Производную в точке всегда легко вычислить, достаточно подставить известные переменные (в момент «настоящее» они известны) в формулу y' = f(x(t), y(t), t). По определению производная связана с углом наклона касательной: y' = tg(α), значит, угол α легко вычислить (α = arctg(y' )) и провести касательную. Проведем касательную до пересечения с линией t + Δt. Момент t + Δt соответствует «будущему» состоянию системы. Проведем линию параллельно оси t от точки A до пересечения с линией t + Δt. Линии образуют прямоугольный треугольник ABC, один катет которого равен Δt (известен). Известен также угол α. Тогда второй катет в прямоугольном треугольнике ABC равен: a = Δt · tg(α). Теперь легко вычислить ординату точки B. Она состоит из двух отрезков — y(t) и a. Ордината символизирует положение системы в точке y(t + Δt). То есть y(t + Δt) = y(t) + a или далее y(t + Δt) = y(t) + Δt · tg(α) или, подставляя дальше: y(t + Δt) = y(t) + Δt · y' и, наконец, y(t + Δt) = y(t) + Δt · f(x(t), y(t), t). Т. е. : будущее = настоящее + изменение 128 будущее = настоящее + скорость · шаг
Метод Эйлера Формула Эйлера может дать точные результаты только при очень малых Δt (при Δt –> 0). При Δt≠ 0 формула дает расхождение между истинным значением y и расчетным, равное ε, поэтому: y(t + Δt) = y(t) + Δt · f(x(t), y(t), t) + ε. Если мысленно сдвигать линию t + Δt на рис. (геометрическая иллюстрация) влево (приближать значение Δt к нулю), расстояние BB* = ε будет сокращаться. Численный метод Эйлера на каждом шаге имеет погрешность расчета ε. Чем меньше взять величину Δt, тем меньше будет ошибка расчета ε (см. рис. ). Для расчета поведения системы на сколько-нибудь продолжительном отрезке времени (например, от t 0 до tk), чтобы уменьшить ошибку на каждом шаге, шаги Δt делают по возможности малыми. Для достижения точки tk отрезок (tk – t 0) делится на отрезки длиной Δt; таким образом, всего получится N = (tk – t 0)/Δt шагов. В результате расчета придется формулу Эйлера применить для каждого шага, то есть N раз. Но следует иметь в виду, что ошибки εi на каждом i-ом шаге (в простейшем случае) складываются, а общая ошибка быстро накапливается. В этом состоит существенный недостаток данного метода. 129
Метод Эйлера Хотя с помощью этого метода можно получить (в численном виде) решение любого дифференциального уравнения (в том числе и неразрешимого аналитически). Уменьшая шаг, мы получаем более точные решения, но при этом не следует забывать, что увеличение числа шагов ведет к вычислительным затратам и снижению быстродействия. Кроме того, при большом числе итераций в расчет вносится другая существенная погрешность из-за ограниченной точности вычислительных машин и ошибок округления. Пример Дано дифференциальное уравнение y' = 2 ty. Задано начальное положение системы: y(0) = 1. Требуется найти y(t), то есть поведение системы на интервале времени t от 0 до 1. Аналитический способ решения y' = 2 ty Методом разделения переменных найдем: y'/y = 2 t Будем интегрировать от 0 до ti, тогда согласно правилам интегрирования имеем: Полученное аналитическое решение характеризуется тем, что оно является абсолютно точным, но если уравнение окажется сколько-нибудь сложным, то решение не будет найдено вовсе. Аналитический путь решения не универсален. 130
Метод Эйлера Пример Численный способ решения предполагает, что расчет будет вестись по формуле Эйлера на ряде последовательных шагов. На каждом шаге решение имеет свою ошибку, поскольку на каждом шаге кривая заменяется прямым отрезком. При алгоритмической реализации расчет реализуется циклом, в котором изменяется t (счетчик t) и y: t : = t + Δt; y : = y + 2 · t · y · Δt Будем искать значение y в численном виде на промежутке от T = 0 до T = 1. Возьмем число шагов n = 10, тогда шаг приращения Δt составит: Δt = (1 – 0)/n = (1 – 0)/10 = 0. 1. Численный расчет уравнения методом Эйлера и сравнение результата с точным решением на каждом шаге i ti 0 1 2 3 4 5 6 7 8 9 10 0. 1 0. 2 0. 3 0. 4 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 yi = yi – 1 + y'i – 1 · Δt 1 1 1. 02 1. 061 1. 124 1. 214 1. 336 1. 496 1. 706 1. 979 2. 335 y'i = 2 ti · yi 0 0. 2 0. 408 0. 636 0. 900 1. 214 1. 603 2. 095 2. 729 3. 561 4. 669 Δyi = y'i · Δt 0 0. 02 0. 0408 0. 0636 0. 0900 0. 1214 0. 1603 0. 2095 0. 2729 0. 3561 0. 4669 yi + 1 = yi + yточн. = exp(ti 2) Δyi 1 1 1. 02 1. 0101 1. 0608 1. 0408 1. 1246 1. 0942 1. 2140 1. 1735 1. 3354 1. 2840 1. 4963 1. 4333 1. 7055 1. 6323 1. 9789 1. 8965 2. 3351 2. 2479 2. 8019 2. 7183 131
Метод Эйлера Пример Численный способ решения Рассчитанное численно значение (yi + 1) отличается от точного (yточн. ), и погрешность (разница столбцов yi + 1 и yточн. ) в процессе расчета нарастает. Относительная погрешность σ для расчетного значения y(1), полученного численно, в сравнении с теоретическим точным yтеор. : σ = (1 – yрасч. /yтеор. ) · 100%. Если менять значение шага Δt, например, уменьшать шаг, то относительная погрешность расчета тоже будет уменьшаться. Δt yрасч. (1) yтеор. (1) σ Зависимость погрешности расчета от размера шага Δt: С уменьшением шага приращения Δt уменьшается 1/10 2. 3346 2. 7183 14% величина относительной погрешности, а значит, повышается точность расчета. 1/20 При этом изменение шага в 10 раз (с 1/10 до 1/100) ведет к изменению величины ошибки примерно тоже в 10 раз (с 14 до 2%). 1/100 2. 5107 2. 6738 2. 7183 8% 2. 7183 2% Размер шага и ошибка для метода Эйлера связаны линейно. Это принято обозначать символом ε = O(Δt), а метод Эйлера называют методом первого порядка точности. Поскольку в методе Эйлера ошибка достаточно велика и от шага к шагу накапливается, а точность пропорциональна количеству вычислений, то метод Эйлера обычно применяют для грубых расчетов, для оценки поведения системы в принципе. Для точных количественных расчетов применяют более точные методы. 132
Метод Эйлера Примечания 1. Каждый численный метод обладает точностью, поскольку результат отличается от теоретического. Точность метода зависит от величины шага. Различные методы имеют различную точность. Порядок зависимости точности от величины шага обозначают как O(h). У метода Эйлера первый порядок точности, зависимость ошибки от величины шага линейна. 2. Если при уменьшении шага предел yn стремится к значению yтеор. , то говорят, что метод обладает сходимостью. Исследователей интересует скорость сходимости метода. 3. Метод должен быть устойчив. Устойчивость связана с некоторой критической величиной шага. При проявлении неустойчивости наблюдается полное искажение качественной картины расчета, «разболтка» результата. 4. При выборе метода рекомендуется сначала добиться устойчивости, а внутри области устойчивости — сходимости результата. Устойчивость обеспечивает качественную картину. Сходимость обеспечивает количественный результат. 133
Метод Эйлера Пример Уравнения описывают процесс теплообмена между двумя телами, температуры которых в некоторый момент времени обозначим как A и B — переменные, меняющиеся во времени t. Найти, как будут меняться температуры A(t) и B(t) (найти поведение системы). Интуитивно ясно, что при начальной разнице температур A = 8 и B = 5 температуры тел постепенно со временем должны выровняться, так как более горячее тело будет отдавать энергию более холодному, и его температура будет уменьшаться, а более холодное тело будет принимать энергию от более горячего, и его температура будет увеличиваться. Процесс теплообмена закончится (то есть изменения прекратятся) тогда, когда температуры двух тел станут одинаковыми. Система ведет себя качественно Проведем несколько расчетов поведения A(t) и B(t) с неверно. Решение неустойчиво разной величиной шага Δt. Будем брать различную величину шага Δt и находить соответствующие значения A и B во времени по следующим формулам Эйлера: Aнов. = Aпред. + (Bпред. – Aпред. ) · Δt, нов. = Bпред. + (Aпред. – Bпред. ) · Δt. B Расчет при Δt = 2: № шага t A B 0 0 8 5 1 2 2 11 2 4 20 – 7 Наблюдается явление «разболтки» . Неустойчивое решение. Из физических соображений очевидно, что так вести себя два тела в процессе теплообмена не могут. 134
Пример Расчет при Δt = 1: Метод Эйлера № шага t A 0 0 8 5 1 1 5 8 2 2 8 Система ведет себя качественно неверно. Решение находится на грани устойчивости B 5 Наблюдается поведение решения системы на границе устойчивости. Расчет при Δt = 0. 5: Система ведет себя качественно правильно. Решение (поведение системы) имеет большую погрешность № шага t A B 0 0 8 5 1 0. 5 6. 5 2 1. 0 6. 5 Решение устойчиво, соответствует правильной качественной картине. Температуры тел постепенно сближаются, становятся со временем одинаковыми. Но решение пока имеет большую погрешность. 135
Пример Расчет при Δt = 0. 1: № шага Решение устойчиво. 0 Решение более точно. 1 2 Метод Эйлера t A B 0 0. 1 0. 2 8 7. 7 7. 46 5 5. 3 5. 54 3 0. 3 7. 27 5. 73 4 0. 4 7. 12 5. 88 5 0. 5 7. 00 Система ведет себя качественно верно. Количественно решение более точно 6. 00 Связь величины шага расчета с устойчивостью метода и его точностью 136
Метод Эйлера Можно ли увеличить точность на порядок, но при этом сэкономить на количестве вычислений? Модифицированный метод Эйлера имеет точность второго порядка. В методе Эйлера производная берется в начале шага и по ней прогнозируется движение системы на конец шага, считая, что во время шага производная неизменна. То есть в течение всего шага производная считается той, какой она была в самом начале шага. Это основной источник неточности. Улучшение метода состоит в том, что берется производная не в начале шага, а как промежуточное или среднее на разных участках одного шага. В разных вариантах метода вычисляют несколько производных в разных частях шага и усредняют их. Конечно, в этом случае число вычислений увеличивается, — но не в десятки раз, — а вот точность возрастает на порядок, в этом и состоит выигрыш. Пример. Пусть (как и прежде) требуется решить уравнение y' = f(x, y, t). Идея уточненного метода Эйлера состоит в том, что производную вычисляют не в i-ой точке, а между двумя соседними точками: i и i + 1. Данная процедура состоит из следующих шагов: 1. в точке i вычисляют значение производной: f(xi, yi, t); 2. делают полшага и вычисляют значение функции на середине отрезка: yi + 1/2 = yi + f(xi, yi, t) · Δt/2; 3. в точке i + 1/2 вычисляют производную: f(xi + 1/2, yi + 1/2, t + Δt/2); делается полный шаг из точки i в точку i + 1 по значению уточненной производной: yi + 1 =yi + f(xi + 1/2, yi + 1/2, t + Δt/2) · Δt; 4. значение t увеличивается: t : = t + Δt. Вся процедура повторяется сначала. Данный метод обладает точностью Ο 2(h), то есть на порядок выше, 137 чем метод Эйлера, при увеличении числа вычислений всего в 2 раза.
Метод Эйлера Какой будет ошибка ε (расхождение между реальным и вычисленным теоретическим значением), если шаг делается по значению производной, вычисленной в точке i, как это делается в методе Эйлера. Эта ошибка может быть достаточно велика! Движение реальной и расчетной системы по методу Эйлера и расхождение между ними (ошибка) 138
Метод Эйлера Модифицированный метод Эйлера По значению производной, вычисленной в точке i, делается полшага до точки t + Δt/2 (направление производной показано линией A). И в точке t + Δt/2 вычисляют новую производную. Касательная в точке t + Δt/2 будет другой — линия B. Ее наклон равен производной в точке t + Δt/2. Уточнение значения производной внутри шага расчета 139
Метод Эйлера Модифицированный метод Эйлера Переносят линию B обратно в точку t. Это соответствует тому, что из точки t снова делается, — но уже полный, — шаг Δt до точки t + Δt по направлению, соответствующему линии C. Линия C параллельна B. То есть значение производной в точке t берется искусственно равным производной в точке t + Δt/2. Ошибка расчета (ε 1) во многих случаях при этом уменьшается. Движение реальной системы и системы, рассчитанной модифицированным методом Эйлера, и расхождение между ними (ошибка). На рисунке для сравнения показан результат, вычисленный по методу Эйлера 140
Метод Эйлера Модифицированный метод Эйлера Существует другой вариант модифицированного метода Эйлера, когда производную для того, чтобы сделать шаг из точки i, берут не в i-ой точке и не в i + 1/2, а как среднее арифметическое двух производных: производной в точке i (направление производной показано линией A) и производной в точке i + 1 (направление производной показано линией B). Направление «средней» производной показано линией C. 141
Метод Рунге-Кутты используют для расчета стандартных моделей достаточно часто, так как при небольшом объеме вычислений он обладает точностью метода Ο 4(h). Для построения разностной схемы интегрирования воспользуемся разложением функции в ряд Тейлора: Заменим вторую производную в этом разложении выражением где Причем Δx подбирается из условия достижения наибольшей точности записанного выражения. Для дальнейших выкладок произведем замену величины «y с тильдой» разложением в ряд Тейлора: Для исходного ДУ построим вычислительную схему: которую преобразуем к виду: 142
Метод Рунге-Кутты Введем следующие обозначения: Эти обозначения позволяют записать предыдущее выражение в форме: Очевидно, что все введенные коэффициенты зависят от величины Δx и могут быть определены через коэффициент α, который в этом случае играет роль параметра: Окончательно схема Рунге-Кутты принимает вид: Та же схема в форме разностного аналога исходного ДУ: При α = 0 получаем как частный случай уже известную схему Эйлера: При α = 1: 143
Метод Рунге-Кутты При α = 1 проведение расчетов на очередном шаге интегрирования можно рассматривать как последовательность нижеследующих операций. 1. Вычисляется выражение, представляющее собой полушаг интегрирования по схеме Эйлера, то есть определяется приближенное значение искомой функции в точке xk + h/2: 2. Для той же промежуточной точки находится приближенное значение производной: 3. Определяется уточненное значение функции в конечной точке всего шага, причем по схеме Эйлера с вычисленным на предыдущем шаге значением производной: Геометрические построения показывают, что получаемое в такой последовательности решение лежит «ближе» к истинному, чем вычисляемое по схеме Эйлера, то есть следует ожидать более высокой точности решения, получаемого методом Рунге-Кутты. Ранее мы назвали эту схему «модифицированным методом Эйлера» . 144
Метод Рунге-Кутты Теперь рассмотрим схему при α = 0. 5: 1. Выполняется полный шаг метода Эйлера с целью определения приближенного значения искомой функции на конце отрезка интегрирования: 2. Для этой же точки вычисляется приближенное значение производной: 3. Находится среднее значение двух производных, определенных на концах отрезка: 4. Вычисляется значение искомой функции в конечной точке всего шага по схеме Эйлера с усредненным значением производной: Иногда получающееся выражение называют схемой (методом) Эйлера-Коши. Геометрически понятно, что получаемый указанным способом результат также должен быть «ближе» к истинному решению, чем получаемый по схеме Эйлера. 145
Метод Рунге-Кутты третьего и четвертого порядков Рассмотрим две различные схемы Рунге-Кутты, предназначенные для численного решения обыкновенных дифференциальных уравнений первого порядка и имеющие третий порядок аппроксимации: И две схемы Рунге-Кутты, имеющие четвертый порядок аппроксимации: Более высокая степень аппроксимации дифференциального уравнения разностным аналогом позволяет получать более точное решение при более крупном шаге и, следовательно, меньшем числе шагов, то есть приводит к снижению требуемых ресурсов ЭВМ. 146 На сегодняшний день для грубого расчета вычисления производятся методом Эйлера, для точного расчета — методом Рунге-Кутты.
Статистическое моделирование — базовый метод моделирования, заключающийся в том, что модель испытывается множеством случайных сигналов с заданной плотностью вероятности. Целью является статистическое определение выходных результатов. В основе статистического моделирования лежит метод Монте-Карло. Имитацию используют тогда, когда другие методы применить невозможно. Метод Монте-Карло Рассмотрим метод Монте-Карло на примере вычисления интеграла, значение которого аналитическим способом найти не удается. Задача 1. Найти значение интеграла: Вычислить значение интеграла функции f(x) – значит найти площадь под графиком. Ограничим кривую сверху, справа и слева. Случайным образом распределим точки в прямоугольнике поиска. Обозначим N 1 количество точек, принятых для испытаний (то есть попавших в прямоугольник, изображены красным и синим цветом), N 2 количество точек под кривой, то есть попавших в закрашенную площадь под функцией (изображены красным цветом). Тогда естественно предположить, что количество точек, попавших под кривую по отношению к общему числу точек пропорционально площади под кривой (величине интеграла) по отношению к площади 147 прямоугольника.
Метод Монте-Карло Значения r 1 и r 2 на являются равномерно распределенными случайными числами из интервалов (x 1; x 2) и (c 1; c 2) соответственно. Метод Монте-Карло чрезвычайно эффективен, прост, но необходим «хороший» генератор случайных чисел. Вторая проблема применения метода заключается в определении объема выборки, то есть количества точек, необходимых для обеспечения решения с заданной точностью. Эксперименты показывают: чтобы увеличить точность в 10 раз, объем выборки нужно увеличить в 100 раз; то есть точность примерно пропорциональна корню квадратному из объема выборки: 148
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Построив модель системы со случайными параметрами, на ее вход подают входные сигналы от генератора случайных чисел (ГСЧ). ГСЧ устроен так, что он выдает равномерно распределенные случайные числа rрр из интервала [0; 1]. Так как одни события могут быть более вероятными, другие — менее вероятными, то равномерно распределенные случайные числа от генератора подают на преобразователь закона случайных чисел (ПЗСЧ), который преобразует их в заданный пользователем закон распределения вероятности, например, в нормальный или экспоненциальный закон. Эти преобразованные случайные числа x подают на вход модели. Модель отрабатывает входной сигнал x по некоторому закону y = φ(x) и получает выходной сигнал y, который также является случайным. В блоке накопления статистики (БНСтат) установлены фильтры и счетчики. Фильтр (некоторое логическое условие) определяет по значению y, реализовалось ли в конкретном опыте некоторое событие (выполнилось условие, f = 1) или нет (условие не выполнилось, f = 0). Если событие реализовалось, то счетчик события увеличивается на единицу. Если событие не реализовалось, то значение счетчика не меняется. Если требуется следить за несколькими разными типами событий, то для статистического моделирования понадобится несколько фильтров и счетчиков Ni. 149 Всегда ведется счетчик количества экспериментов — N.
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Далее отношение Ni к N, рассчитываемое в блоке вычисления статистических характеристик (БВСХ) по методу Монте-Карло, дает оценку вероятности pi появления события i, то есть указывает на частоту его выпадения в серии из N опытов. Это позволяет сделать выводы о статистических свойствах моделируемого объекта. Например, событие A совершилось в результате проведенных 200 экспериментов 50 раз. Это означает, согласно методу Монте-Карло, что вероятность совершения события равна: p. A = 50/200 = 0. 25. Вероятность того, что событие не совершится, равна, соответственно, 1– 0. 25 = 0. 75. Обратите внимание: когда говорят о вероятности, полученной экспериментально, то ее называют частостью; слово вероятность употребляют, когда хотят подчеркнуть, что речь идет о теоретическом понятии. При большом количестве опытов N частота появления события, полученная экспериментальным путем, стремится к значению теоретической вероятности появления события. В блоке оценки достоверности (БОД) анализируют степень достоверности статистических экспериментальных данных, снятых с модели (принимая во внимание точность результата ε, заданную пользователем) и определяют необходимое для этого количество статистических испытаний. Если колебания значений частоты появления событий относительно теоретической вероятности меньше заданной точности, то экспериментальную частоту принимают в качестве ответа, иначе генерацию случайных входных воздействий продолжают, и процесс моделирования повторяется. При малом числе испытаний результат может оказаться недостоверным. Но чем более испытаний, тем точнее ответ, согласно центральной предельной теореме. Заметим, что оценивание ведут по худшей из частот. Это обеспечивает достоверный результат сразу по всем снимаемым характеристикам модели. 150
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Пример. Какова вероятность выпадения монеты орлом кверху при падении ее с высоты случайным образом? Начнем подбрасывать монетку и фиксировать результаты каждого броска. Количество опытов N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Значение счетчика выпадения орла Nо 0 0 1 1 2 3 4 … … … … Значение счетчика выпадения решки Nр 1 2 2 3 3 … … … … Частость выпадения орла Pо =Nо/N 0 0 0. 33 0. 25 0. 4 0. 57 … … … … Частость выпадения решки Pр =Nр/N 1 1 0. 66 0. 75 0. 6 0. 5 0. 43 … … … … Будем подсчитывать частость выпадения орла как отношение количества случаев выпадения орла к общему числу наблюдений. Случаи для N = 1, N = 2, N = 3 — сначала значения частости нельзя назвать достоверными. 151
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами График зависимости Pо от N: как меняется частость выпадения орла в зависимости от количества проведенных опытов (при различных экспериментах будут получаться разные таблицы и, следовательно, разные графики). 152
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Некоторые выводы: 1. Видно, что при малых значениях N, например, N = 1, N = 2, N = 3 ответу вообще доверять нельзя. Например, Pо = 0 при N = 1, вероятность выпадения орла при одном броске равна нулю! Хотя всем хорошо известно, что это не так. То есть пока мы получили очень грубый ответ. Однако в процессе накопления информации ответ медленно, но верно приближается к правильному. К счастью, в данном случае правильный ответ известен: в идеале, вероятность выпадения орла равна 0. 5 (в других более сложных задачах ответ, конечно, будет неизвестен). Допустим, что ответ надо знать с точностью ε = 0. 1. Проведем две параллельные линии, отстоящие от правильного ответа 0. 5 на расстояние 0. 1. Ширина образовавшегося коридора будет равна 0. 2. Как только кривая Pо(N) войдет в этот коридор так, что уже его не покинет, можно остановиться и посмотреть, для какого значения N это произошло. Это и есть экспериментально вычисленное критическое значение необходимого количества опытов Nкрэ для определения ответа с точностью ε = 0. 1; ε-окрестность в этих рассуждениях играет роль своеобразной трубки точности. Заметьте, что ответы Pо(91), Pо(92) и так далее уже не меняют сильно своих значений; по крайней мере, у них не изменяется первая цифра после запятой, которой мы обязаны доверять по условиям задачи. 2. Причиной такого поведения кривой является действие центральной предельной теоремы. В самом простом варианте «Сумма случайных величин есть величина неслучайная» . Использована средняя величина Pо, которая несет в себе информацию о сумме опытов, и поэтому постепенно эта величина становится все более достоверной. 3. Если проделать еще раз этот опыт сначала, то, конечно, его результатом будет другой вид случайной кривой. И ответ будет другим, хотя примерно таким же. Проведем серию таких экспериментов. Такая серия называется ансамблем реализаций. Какому же ответу в итоге следует верить? Ведь они, хоть и являются близкими, все же разнятся. На практике поступают поразному. Первый вариант — вычислить среднее значение ответов за несколько реализаций. 153
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами 154
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Экспериментальные данные необходимого количества бросков монеты для достижения точности ε = 0. 1 при вычислении вероятности выпадения орла Опыт Nкрэ 1 288 2 95 3 50 4 29 5 113 6 210 7 30 8 42 9 39 10 48 Среднее Nкр. э 94 Мы поставили несколько экспериментов и определяли каждый раз, сколько необходимо было сделать опытов, то есть Nкрэ. Было проделано 10 экспериментов, результаты которых были сведены в табл. По результатам 10 -ти экспериментов было вычислено среднее значение Nкрэ. Таким образом, проведя 10 реализаций разной длины, мы определили, что достаточно в среднем было сделать 1 реализацию длиной в 94 броска монеты. На рис. представлено 100 реализаций — 100 красных линий. Абсцисса N = 94 отмечена вертикальной чертой. Есть какой-то процент красных линий, которые не успели пересечь εокрестность, то есть (Pэксп – ε ≤ Pтеор ≤ Pэксп + ε), и войти в коридор точности до момента N = 94. Таких линий 5. Это значит, что 95 из 100, то есть 95%, линий достоверно вошли в обозначенный интервал. Таким образом, проведя 100 реализаций, мы добились примерно 95%-ного доверия к полученной экспериментально величине вероятности выпадения орла, определив ее с точностью 0. 1. 155
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Теоретический расчет необходимого количества бросков монеты для достижения точности ε = 0. 1 при вычислении вероятности выпадения орла Доверительная вероятность QF 0. 90 0. 95 0. 99 Коэффициент Лапласа k(QF) 2. 72 3. 84 6. 66 Требуемое число опытов Nкрт = k(QF) · p · (1 – p)/ε 2 68 96 167 Для сравнения полученного результата вычислим теоретическое значение Nкрт теоретически. Однако для этого придется ввести понятие доверительной вероятности QF, которая показывает, насколько мы готовы верить ответу. Например, при QF = 0. 95 мы готовы верить ответу в 95 % случаев из 100. Формула теоретического расчета числа экспериментов имеет вид: Nкрт = k(QF) · p · (1 – p)/ε 2, где k(QF) — коэффициент Лапласа, p — вероятность выпадения орла, ε — точность (доверительный интервал). В табл. показаны значения теоретической величины количества необходимых опытов при разных QF (для точности ε = 0. 1 и вероятности p = 0. 5). Как видите, полученная нами оценка длины реализации, равная 94 опытам очень близка к теоретической, равной 96. Некоторое несовпадение объясняется тем, что, видимо, 10 реализаций недостаточно для точного вычисления Nкрэ. Если вы решите, что вам нужен результат, которому следует доверять больше, то измените значение доверительной вероятности. Например, теория говорит нам, что если опытов будет 167, то всего 1 -2 линии из ансамбля не войдут в предложенную трубку точности. Но имейте в виду, количество экспериментов с ростом точности и достоверности растет очень быстро. 156
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Второй вариант, используемый на практике — провести одну реализацию и увеличить полученное для нее Nкрэ в 2 раза. Это считают хорошей гарантией точности ответа. 157
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Если присмотреться к ансамблю случайных реализаций, то можно обнаружить, что сходимость частости к значению теоретической вероятности происходит по кривой, соответствующей обратной квадратичной зависимости от числа экспериментов. 158
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами Это действительно так получается и теоретически. Если изменять задаваемую точность ε и исследовать количество экспериментов, требуемых для обеспечения каждой из них, то получится табл. График зависимости Nкрт(ε) Теоретическая зависимость количества экспериментов, необходимых для обеспечения заданной точности при QF = 0. 95 Точность ε 0. 1 0. 001 Критическое число экспериментов Nкрт 96 960000 159
Экстраполяция результатов моделирования с животных на человека В процессе исследования влияния факторов среды на человека возможны серьезные ограничения в проведении экспериментов и наблюдений. Эти ограничения обусловлены прежде всего требованиями техники безопасности. Наряду с этим могут возникнуть трудности технического и материального характера. В таких случаях эксперименты целесообразно проводить на животных с последующим переносом результатов анализа полученных данных на человека. Задачи • Выбор среди различных видов животных приближенных биологических моделей человека в отношении реакции на исследуемые факторы среды; • Проведение экспериментов с животными и построение математических моделей для каждой биологической модели; • Построение масштабов, связывающих сходственные параметры математических моделей для каждой пары биологических моделей; • Построение функциональных зависимостей, связывающих величины масштабов для каждого сходственного параметра с характеристиками организма; • Экстраполяция результатов моделирования в ряду биологических моделей 160 на человека.
Выбор биологических моделей Биологические модели человека – это млекопитающие различных видов, имеющие общность фундаментальных физико-химических процессов, протекающих в организме, по отношению к исследуемому влиянию факторов среды. Эти процессы в связи с их большой сложностью математически описываются обычно лишь приближенно. Поэтому необходимое условие выбора животных в качестве биологических моделей состоит в том, чтобы их реакции на комплекс изучаемых воздействий были качественно сходными с реакциями человека. В связи с этим животные могут рассматриваться только как приближенные биологические модели человека. Например, белые мыши являются биологическими моделями человека при исследовании и оптимизации лечения ряда инфекционных заболеваний ( грипп, венесуэльский энцефалит и т. п. ) Для решения задачи экстраполяции результатов моделирования на человека целесообразно выбрать несколько биологических моделей ( не менее трех). Такое условие обусловлено необходимостью построения функциональных зависимостей, связывающих между собой сходственные параметры математических моделей реакций выбранных биологических объектов. 161
Проведение экспериментов с животными и построение математических моделей Вид эксперимента и структура математической модели зависят от цели исследования, количества воздействующих факторов среды и методов оценки состояния организма. Так, в случае оценки реакции организма по отдельным показателям в ответ на комбинированное воздействие факторов среды целесообразно использовать модели «поведения» , построенные на основе методологии планирования эксперимента: где Yk — функция отклика k-го вида животных; bi, 0, ii — коэффициенты модели; n — количество воздействующих факторов. Если реакция организма оценивается по альтернативной классификации, например на основе методов экспертных оценок, то наиболее адекватно построение моделей на основе теории распознавания образов: где Fk(A, X) — уравнение гиперплоскости для k-го вида животных; a 0 i — коэффициенты модели; n – количество коэффициентов. В результате выполнения данного этапа экспериментального исследования получается система из k экспериментально-статистических моделей, каждая из которых должна быть адекватна соответствующей биологической модели. Для продолжения исследований необходимо, чтобы знаки перед наиболее значащими коэффициентами совпадали, так как в противном случае не выполняется условие качественного сходства моделей. 162
Построение масштабов, связывающих сходственные параметры математических моделей для каждой пары биологических моделей Для реализации экстраполяции или переноса результатов математического моделирования с выбранных биологических моделей на человека необходимо осуществить процедуру, связывающую между собой сходственные параметры построенных моделей. Эта процедура основана на некоторых элементах теории подобия, устанавливающих функциональные связи между параметрами процессов, протекающих в оригинале и модели, расчетом масштабов или констант подобия. Сходственные параметры математических моделей – это переменные, имеющие одинаковые индексы и соответствующие подобным характеристикам объектов или воздействующим факторам. Подобие – это взаимооднозначное соответствие между двумя объектами, при котором функции перехода от параметров, характеризующих один из объектов, к сходственным параметрам другого объекта известны, а математические описания этих объектов могут быть преобразованы в тождественные. Константа подобия (М) – это некоторая величина, постоянная для каждой пары биологических объектов, один из которых рассматривается как оригинал, а другая – как модель: Mj = Xj /xj, j = 1, 2, …, k, где Xj – j-й параметр оригинала; xj – j-й параметр модели; k – количество параметров модели. В качестве параметров моделей могут рассматриваться воздействующие факторы, показатели реакции объектов на воздействия, характеристики структуры объектов ( для биологии – анатомические, физиологические, морфологические, биохимические). При выполнении условия подобия становится возможным, исследовав реакции модели на действие комплекса факторов, оценить реакции оригинала на действие 163 подобных факторов: Xj = Mjxj, j = 1, 2, …, k.
Построение масштабов, связывающих сходственные параметры математических моделей для каждой пары биологических моделей Например для двух видов животных, в эксперименте получены две математические модели, связывающие показатели состояния Y и y с комбинациями воздействующих факторов, соответственно, (X 1, X 2, …, Xn ) для Y и (x 1, x 2, …, xn) для у: Y=F(X 1, X 2, …, Xn), y=f(x 1, x 2, …, xn), где n – количество воздействующих факторов. Первая зависимость описывает реакцию оригинала (Y), а вторая – реакцию модели (у). Необходимо определить масштабы: My=Y/y, Mj=Xj /xj, j=1, 2, …, n. Из предложения подобия следует, что существует такая комбинация масштабов My* и Mj*, которая обращает приведенные уравнения в тождество: Y=y. При этом функции F и f в обоих уравнениях предполагаются одинаковыми. Для определения масштабов целесообразно реализовать процедуру оптимизации с целевой функцией С, характеризующей отклонения значений реакции у модели, умноженной на масштаб My, от реакции оригинала Y. Реакции моделей вычисляются в заданных точках пространства х1, х2, …, хn, определяемых, например планами 2 -го порядка, в соответствии с требованиями методологии факторного эксперимента. Реакции оригинала рассчитываются в сходственных точках этого же пространства, координаты которого равны Xj =Mjxj, где j=1, 2, …, n. Среди различных целевых функций наиболее эффективно может быть использована сумма квадратов отклонений реакций оригинала и модели в сходственных точках: где N – количество комбинаций величин факторов х1, х2, …, хn, задаваемых планом эксперимента. 164
Построение масштабов, связывающих сходственные параметры математических моделей для каждой пары биологических моделей Масштабы My и Mi, i=1, 2, …, n, являются искомыми величинами, минимизирующими целевую функцию. Для отыскания масштабов может быть использована любая процедура поисковой оптимизации. Аналогично рассчитываются масштабы для многомерных моделей альтернативной классификации в виде уравнений гиперплоскости. В таком случае строятся два уравнения, соответственно, для объекта-оригинала (f 1) и его модели (f 2) вида где n – количество сходственных параметров (Xi и xi) могут рассматриваться показатели состояния исследуемых биологических объектов и воздействующие факторы, например физическая нагрузка или доза фармакологического средства. В процессе распознавания, т. е. отнесения объекта к одной из двух групп, необходимо вычислить величину и знак отклонения соответствующей ему точки гиперплоскости в пространстве сходственных параметров ( Y – для оригинала и y – для модели): Далее процедура осуществляется так же, как для модели «поведения» в виде полинома. 165
Построение функциональных зависимостей, связывающих масштабы сходственных параметров с размерами тела Из сравнительной физиологии известно, что многие параметры организма животных различных видов (х) связаны с размерами тела (q) зависимостью x=aqb, где а и b – постоянные величины. В качестве показателей размера тела могут рассматриваться, например, масса, площадь поверхности, масса мозга и т. д. поэтому можно связать масштаб каждого сходственного параметра с выбранным показателем размера тела зависимостью или какой-либо другой функцией одной переменной: Mi =f(q), где Mi – масштаб i-го сходственного параметра. При этом уравнение легко линеаризуется логарифмированием обеих его частей и приводится к следующему виду: Для построения таких зависимостей необходимо провести эксперименты, по крайней мере, на трех видах животных с различными значениями q. 166
Экстраполяция математической модели реакции с биологических моделей на человека В результате постановки соответствующего показателя размера тела человека (q) в функциональную зависимость можно вычислить масштаб i-го сходственного параметра человека (Mi) по отношению к одному из исследованных видов животных, который рассматривается как базовый. С помощью рассчитанных масштабов становится возможным рассчитать сходственные параметры и построить модель реакции организма человека. Для полиномиальных моделей Модель реакции организма человека: Модель реакции организма животного: Если принять условие точного подобия этих моделей, то после деления на My Каждого слагаемого модели реакции животного можно получить модель реакции организма человека в окончательном виде: Условие точного подобия представляется в виде констант B 0/Mybo=1, Bi. Mi/Mybi=1, Bij. Mi. Mj/Mybij=1, Bii. Mi 2 /Mybii=1, где B 0, Bij, Bii – коэффициенты модели реакции животного; bo, bij, bii – коэффициенты модели реакции человека; Mi, Mj – масштабы. Т. о. любой коэффициент модели реакции организма человека можно рассчитать 167 следующим образом: bo=B 0/My, bi=Bi. Mi/My, bij=Bij. Mi. Mj/My, b=Bii. Mi 2/My.
Экстраполяция математической модели реакции с биологических моделей на человека Пример: необходимо исследовать эффективность нового профилактического препарата Z против особо опасной инфекции W. Очевидно, что эксперименты с человеком в данном случае невозможны. Поэтому была поставлена задача провести многофакторные эксперименты на трех видах животных, у которых основные качественные показатели развития исследуемой инфекции сходны. По результатам этих экспериментов должны быть построены полиномиальные модели реакций, рассчитаны масштабы и модель реакции для человека. Были выбраны следующие биологические виды: белые мыши, морские свинки и собаки. На каждом из этих видов животных поставлены двухфакторные эксперименты 1 -го порядка. В качестве воздействующих факторов рассматривались доза возбудителя (х1) и доза препарата (х2). Оценивался процент выживаемости животных (у). Условное обозначение Наименование Х 1 Х 2 Диапазоны варьирования величин воздействующих факторов Значение факторов Единица измерения min max Доза возбудителя 0, 2 1 LD 50 Доза препарата 0, 2 1 ED 50 LD 50 – доза возбудителя, при которой погибают 50 % экспериментальных животных; ED 50 – доза препарата, при которой выживают 50 % животных на фоне заражения дозой 2 LD 50 168
Скачать презентацию Санкт-Петербургский Государственный Электротехнический Университет ЛЭТИ Кафедра БТС Моделирование
Моделирование.ppt