Типовые математические схемы моделирования. Примеры и задачи

Л. П. МОХРАЧЕВА ТИПОВЫЕ МАТЕМАТИЧЕСКИЕ СХЕМЫ МОДЕЛИРОВАНИЯ. ПРИМЕРЫ И ЗАДАЧИ Учебное пособие Министерство образования и науки Российской Федерации Уральский федеральный университет имени первого Президента России Б. Н. Ельцина Л. П. Мохрачева ТИПОВЫЕ МАТЕМАТИЧЕСКИЕ СХЕМЫ МОДЕЛИРОВАНИЯ. ПРИМЕРЫ И ЗАДАЧИ Учебное пособие Рекомендовано методическим советом Уральского федерального университета для студентов направления подготовки 08.04.01 — Строительство Екатеринбург Издательство Уральского университета 2018 УДК 510.67(075.8) ББК 22.12я73 М86 Рецензенты: М. П. Кащенко, д‑р физ.-мат. наук, проф., завкафедрой физики Уральского государственного лесотехнического университета; А. П. Танкеев, д‑р физ.-мат. наук, проф., главный научный сотруд‑ ник ИФМ УрО РАН Научный редактор — С. И. Тарлинский, канд. физ.-мат. наук, доц. кафедры высшей математики ИнФО УрФУ Мохрачева, Л. П. М86     Типовые математические схемы моделирования. Примеры и за‑ дачи : учебное пособие / Л. П. Мохрачева. — Екатеринбург : Изд-во Урал. ун-та, 2018. — 144 с. ISBN 978-5-7996-2362-3 Рассмотрены математические схемы некоторых детерминированных и стохастических математических моделей. Для каждой из схем приведе‑ ны примеры математических моделей из различных предметных областей и их решения. Использованы аналитические и численные методы получе‑ ния результатов. Представлены компьютерные программы для численных расчетов и задачи для самоконтроля. Библиогр.: 13 назв. Рис. 21. Табл. 14. Прил. 1. УДК 510.67(075.8) ББК 22.12я73 ISBN 978-5-7996-2362-3 © Уральский федеральный университет, 2018 Предисловие Н астоящее пособие предназначено для первоначально‑ го знакомства с математическими методами модели‑ рования. Объем пособия соответствует семестрово‑ му курсу «Методы математического моделирования» (18 часов лекций и 18 часов практических занятий), который читается автором магистрантам строительного института, обучающим‑ ся по различным образовательным программам. Ограниченный объем курса стал причиной выбора только нескольких типовых математических схем, знакомство с ко‑ торыми может представлять практическую ценность для маги‑ странтов. Ориентация на практическое использование знаний и навыков, которые могут быть получены при изучении кур‑ са, определяет структуру изложения материала в пособии. Из‑ учение типовых моделей происходит в следующей последова‑ тельности: 1) знакомство с математическим аппаратом, использован‑ ным для описания модели; 2) изучение методов решения математической модели; 3) программная реализация алгоритмов решения в рамках доступных математических пакетов; 4) анализ полученных результатов и планирование вычис‑ лительного эксперимента. 3 Предисловие В пособии представлены примеры конкретных расчетов для непрерывных и дискретных детерминированных моделей, неко‑ торых моделей оптимизации, непрерывных стохастических мо‑ делей, а также статистических регрессионных моделей. В конце каждого раздела приведены задачи для самокон‑ троля. Автор надеется, что пособие будет полезным студентам и ма‑ гистрантам при изучении и практическом применении мето‑ дов математического моделирования в различных предметных областях. 4 1. Понятие математической модели. Некоторые типовые схемы математического моделирования 1.1. Формальная математическая модель объекта моделирования Рассмотрим схему Y(y1, y2, ..., yn) V(v1, v2, ..., vl) H(h1, h2, ..., hm) X(x1, x2, ..., xk), где прямоугольником обозначен объект (процесс) моделирова‑ ния, X ( x1 , x2 , ..., xk ) — множество внешних воздействий на мо‑ делируемый объект, V (v1 , v2 , ..., vl ) — множество воздействий внешней среды, H ( h1 , h2 , ..., hm ) — множество внутренних со‑ стояний объекта, Y ( y1 , y2 , ..., yn ) — множество выходных пара‑ метров (реакций) моделируемого объекта. Обычно множество внешних воздействий X ( x1 , x2 , ..., xk ) содержит управляемые па‑ раметры. Однако в общем случае все множества Х, Y, V, H могут содержать как детерминированные, так и стохастические пере‑ менные. Рассматривая    эти  множества как некоторые   формаль‑   ные векторы X , V , H , Y , запишем равенствоY = F ( X , V , H , t ), где t — время, F — некоторый оператор, определение структуры которого и является целью моделирования. Приведенное выше равенство является формальной математической моделью. Если 5 1. Понятие математической модели. Некоторые типовые схемы математического моделирования в нем отсутствует время, то модель называется статической, в от‑ личие от приведенной выше динамической модели. Структура оператора F может быть аналитической, численной или стати‑ стической. В связи с этим математические модели можно разде‑ лить на аналитические, численные и статистические. 1.2. Общая схема создания математической модели и работы с ней Математическая модель создается на основе содержательной модели, которая отражает существенные свойства моделируе‑ мого объекта или процесса и содержит описание входных па‑ раметров (управляемых или стохастических), внутренних па‑ раметров и выходных параметров. Создание математической модели происходит путем формализации содержательной моде‑ ли в виде уравнений, систем уравнений различных видов, нера‑ венств или отношений. На этом этапе можно использовать ти‑ повые математические схемы моделирования. Дальнейшая работа с созданной моделью содержит: 1) ее качественный анализ с целью разработки методов ре‑ шения или возможных упрощений, а также тестирования для частных случаев с использованием точных решений; 2) обоснование выбора конкретного алгоритма для получе‑ ния количественных результатов, в том числе обязательно оцен‑ ку точности вычислительного метода. Реализация выбранных численных методов требует создания компьютерных программ. Для не слишком сложных моделей тексты этих программ можно получить с использованием ма‑ тематических пакетов, таких как MathCAD или Mathematica. Наличие конкретных компьютерных программ делает воз‑ можным проведение вычислительного эксперимента, целью ко‑ торого является проверка адекватности данной модели и полу‑ чение дополнительной информации. 6 1.3. Типовые математические схемы моделирования 1.3. Типовые математические схемы моделирования В пособии рассматриваются следующие типовые математи‑ ческие схемы моделирования: 1) непрерывно-детерминированные D‑схемы. Эти схемы применяются для описания различных моделей в теории управления. Математической моделью является задача Коши для обыкновенного дифференциального уравнения или системы обыкновенных дифференциальных уравнений, а также уравне‑ ний в частных производных с различными начальными и крае‑ выми условиями. В моделях оптимального управления и неко‑ торых других моделях оптимизации решения находятся при ограничениях на функции, которые имеют линейный, нели‑ нейный или дифференциальный характер; 2) дискретно-детерминированные F‑схемы (конечные авто‑ маты). С помощью этих схем описываются модели устройств контроля и управления, имеющие дискретный характер работы во времени. Математическая модель состоит из задания началь‑ ного состояния автомата и из уравнений, задающих значения выходных параметров в данный момент времени в зависимо‑ сти от значений входных параметров и внутренних состояний в данный или предшествующий моменты. Вместо уравнений могут быть использованы таблицы или графы; 3) стохастические модели. В стохастических моделях все или часть переменных множеств Х, Y, Н являются случайными ве‑ личинами. В пособии рассматриваются модели регрессионного типа, имеющие большое значение при статистическом анализе наблюдений и непрерывно-стохастические модели (Q‑схемы), которые применяются для описания систем массового обслу‑ живания. 7 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы) 2.1. Эволюционные модели П римерами использования этих схем могут служить мо‑ дели эволюции биологических систем, статические и динамические модели математической физики, а также модели процессов в конкретных предметных областях и некоторые модели оптимизации (например, оптимального управления). Модели оптимизации будут рассмотрены в сле‑ дующем разделе. ПРИМЕР 1. Простейшей эволюционной моделью является модель динамики популяции (модель Мальтуса). Математическая модель Математической моделью является задача Коши для линей‑ ного дифференциального уравнения первого порядка dN (t ) = ( a (t ) - b (t )) N (t ) , N ( 0 ) = N 0 , dt где N (t ) — численность популяции, a (t ) — коэффициент раз‑ множения, b (t ) — коэффициент смертности. Данное уравнение допускает разделение переменных. 8 2.1. Эволюционные модели Решение Разделяем переменные dN (t ) = ( a (t ) - b (t )) dt . Интегрируя N (t ) правую и левую части уравнения, получаем общий интеграл dN (t ) т N (t ) = ln N (t ) = т ( a (t ) - b (t )) dt + ln C , C № 0. Таким образом, общее решение имеет вид N (t ) = Ce т (a (t ) -b(t ))dt . Из начального условия следует N 0 = Ce т (a (t ) -b(t ))dt t =0 - т ( a (t ) - b(t ))dt t =0 . и C = N 0e Если коэффициенты a (t ) и b (t ) являются константами, то ре‑ шение имеет вид N (t ) = N 0e ( a - b )t . Анализ результата Если a > b, то численность популяции неограниченно возрас‑ тает. При a < b численность снижается, и при a = b численность остается постоянной. Если коэффициенты a (t ) и b (t ) зависят от времени, то картина может быть существенно иной. Напри‑ ж pц мер, если a (t ) = cos t , b (t ) = sin з t + ч , то решение имеет вид и 4ш N (t ) = N 0e ж pц жpц sin t - cos з t + ч + cos з ч и 4ш и4ш . На рис. 2.1 представлен характер этой зависимости при N 0 = 50. 9 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы) 3 1ґ10 800 N ( t) 600 400 200 20 40 60 t Рис. 2.1 Ограничения модели Модель не учитывает ограниченность ресурсов и факторы конкуренции с другими популяциями, хотя частично эти фак‑ торы можно учесть, выбирая функции a (t ) и b (t ). Более точная модель была предложена Ферхюльстом (см. задачу № 1 в кон‑ це раздела). ПРИМЕР 2. Модель конкуренции популяций хищников и жертв (модель Лотки — Вольтерра) [3]. Математическая модель Математической моделью является система нелинейных дифференциальных уравнений dN (t ) dt dM (t ) dt = ( a (t ) - b (t ) M (t )) N (t ) , = ( - g (t ) + d (t ) N (t )) M (t ) , N ( 0 ) = N 0, M ( 0 ) = M 0. 10 (2.1) 2.1. Эволюционные модели N (t ) і 0, M (t ) і 0 — число жертв и хищников, соответствен‑ но, a (t ) і 0, b (t ) і 0 — коэффициенты рождаемости и убыва‑ ния жертв за счет хищников, g (t ) і 0, d (t ) і 0 — коэффициенты убывания и роста хищников. Далее будем считать все коэффи‑ циенты постоянными. Анализ модели Из системы следует dN (t ) dM (t ) = (a - bM (t )) N (t ) ( - g + dN (t )) M (t ) или ( - g + dN (t )) dN (t ) = (a - bM (t )) dM (t ) . N (t ) M (t ) Интегрируя обе части, получим - g ln N (t ) + dN (t ) = a ln M (t ) - bM (t ) + C или ( ) ln M a (t ) N g (t ) = dN (t ) + bM (t ) + ln C , C № 0 . Таким образом, имеем первый интеграл в виде M a (t ) N g (t ) = Ce dN (t ) + bM (t ) , C > 0. Из начальных условий получаем M 9a N 0g = Ce dN 0 +bM 0 , т. е. C = M 9a N 0g e - dN 0 -bM 0 . Система (2.1) имеет стационарные решения. Если коэффи‑ dN (t ) dM (t ) = 0, = 0 следует, циенты постоянны, то из условий dt dt g что равновесные значения хищников и жертв равны N = d a иM = . b 11 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы) Решение Так как система (2.1) нелинейная, то N (t ) и M (t ) находим, решая систему численно с использованием процедуры rkfixed (см. Прил.). В приведенном ниже листинге программы исполь‑ зованы обозначения N (t ) = N 0 (t ) и M (t ) = N 1 (t ). Листинг программы. a := 2 b := 0.02 g := 3 d := 0.05 N 0 := 80 N 1 := 30 t max := 10 i := 100 ж a Ч N 0 - b Ч N 0 N1 ц жN0 ц D (t , N ) := з ч R := rkfixed[з ч , 0, t max, D ] и N1 ш и - g Ч N1 + d Ч N 0 N1 ш t := R N := R 1 M := R . 2 Графики полученных решений представлены на рис. 2.2. На рис. 2.3 приведена фазовая траектория. N M 300 300 200 200 M 100 100 2 6 4 t Рис. 2.2 8 10 50 100 150 N Рис. 2.3 Анализ результатов При выбранных значениях коэффициентов равновесные значения равны N = 60, M = 100. Решения имеют периодиче‑ ский характер. Колебания происходят относительно равновес‑ ных значений. Аналогичная модель используется и для описания других про‑ цессов соперничества (см. задачи № 2 и 3 в конце данного раздела). 12 2.2. Некоторые модели математической физики 2.2. Некоторые модели математической физики Рассмотрим простейшую гидродинамическую модель очист‑ ки воды методом осаждения. В основе модели лежит описание движения твердой частицы в жидкости. ПРИМЕР 3. Математическая модель Дифференциальное уравнение движения твердой частицы в жидкости получают из баланса действующих на нее сил тяже‑ сти, трения и подъемной силы. Уравнение движения имеет вид U 2 (t ) dU W - (rТ + xr ж )V = 0, (2.2) 2 dt где r ж — плотность жидкости; V — объем частицы; V — коэф‑ фициент гидравлического сопротивления; x — коэффициент присоединенной массы; W — миделево сечение; U (t ) — ско‑ рость частицы. Для этого уравнения ставится задача Коши с U (t0 ) = U 0 . g (rТ - r ж )V - Vr ж Анализ модели Запишем уравнение в виде dU (t ) A + BU 2 (t ) = C , dt Vr W где A = (rт + xr ж )V > 0, B = ж > 0, С = g (rт - r ж )V > 0. 2 Разделив на А, получим dU (t ) dt + bU 2 (t ) = c, b = B C ,c = . A A 13 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы) Если коэффициенты b и c не зависят от времени, то это урав‑ нение с разделяющимися переменными, в противном случае это специальное уравнение Риккати. Решение Считая b и c постоянными, разделяем переменные или dU c -U 2 b dU = dt c - bU 2 = bdt . Интегрируя, получим с b = bt + b ln C . c c Ub Таким образом, имеем 1 b ln 2 c U+ c b = Ce 2 cbt . c Ub Из начальных условий U ( 0 ) = U 0 следует U+ c b C= . c U0 b Окончательно для решения получаем выражение U 0+ cж c ц 2 cbt c +U 0 ззU 0 + чe bи b чш b U (t ) = . c ж c ц 2 cbt U0 - зU 0 + чe b зи b чш 14 2.2. Некоторые модели математической физики Полученное выражение можно преобразовать, если разде‑ лить числитель и знаменатель сначала на e cbt и затем на e cbt + e - cbt . В результате получим U (t ) = ( cbt ) . cb + bU th ( cbt ) U 0 cb + c Ч th Анализ решения th Так как tlim ®+Ґ ( ) cbt = lim t ®+Ґ e cbt - e- cbt e cbt + e- cbt = 1, то скорость частицы стремится к постоянному значению U (t ) = U 0 cb + c cb + bU 0 . Конкретные значения параметров, входящих в уравне‑ жмц ж кг ц ние (2.2), имеют размерности g з 2 ч, rТ , r ж з 3 ч, V, x — безраз‑ ис ш им ш V мерные величины. Значения варьируются в пределах: 1) тонкая пластинка V » 1,1; 2) открытая полусфера V » 0, 35; 3) сфера V » 0,1 – 0, 4; тело обтекаемой формы V » 0, 04 . Для сфе‑ ры x » 0, 5. Размерности W и V м 2 и м 3, соответственно. Если ча‑ 4 стица имеет форму эллипсоида вращения, то W = pa 2 ,V = pa 2c. 3 На рис. 2.4 представлен график решения U (t ), полученного при следующих значениях параметров:V = 0, 0000001, W = 0, 00003, V = 0, 4, x = 0, 5, rТ = 2400, r ж = 998, 2, g = 9, 8. Анализ модели Модель не учитывает некоторые важные эффекты, которые возникают при осаждении, в частности, возможное увеличе‑ ние массы частицы из-за присоединения других частиц, нали‑ 15 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы) чие в жидкости загрязнений, а также возникновение градиента концентрации, что изменяет условия осаждения. Математиче‑ ская модель, учитывающая эти эффекты, имеет форму общего уравнения Риккати с переменными коэффициентами, которое допускает только численное решение (см. задачу № 4 в конце этого раздела). 0.45 U( t) 0.4 0.35 0.3 0.2 0.4 0.6 0.8 t Рис. 2.4 ПРИМЕР 4. Аналогичные модели в D‑схеме. Эти математические модели описывают процессы различной физической природы при том, что математические модели ана‑ логичны. Этот факт является следствием того, что различные физические законы, лежащие в основе математических моде‑ лей, имеют одинаковое математическое описание: 1. Закон Фурье распространения тепла  Q = -lgradT ,  где Q — поток тепла; Т — температура; λ — коэффициент те‑ плопроводности. 16 2.2. Некоторые модели математической физики 2. Закон Ома  1 J = - gradU , r  где J — поток тока; U — потенциал; ρ — коэффициент электро‑ проводности. 3. Закон Дарси — Вейсбаха  V = -kgradP ,  гдеV — вектор скорости; P — давление; k — коэффициент филь‑ трации. 4. Закон Фика  M = -DgradC ,  где M – поток массы; С — концентрация; D — коэффициент диффузии. В качестве примера рассмотрим уравнение теплопроводности ¶U ( M , t ) = div ( lgradU ) + F ( M ,t ) , ¶t где M ( x, y, z ) — точка в пространстве; U ( M , t ) — температура; cr с — удельная теплоемкость; ρ — плотность вещества; λ — коэф‑ фициент теплопроводности; F ( M , t ) — плотность источников тепла. Далее будем считать, что с, ρ и λ — постоянные, и преобра‑ зуем уравнение к виду ¶U ( M , t ) ¶t где a = = aDU ( M , t ) + 1 F (M , t ), cr ¶2 ¶2 ¶2 l ; D = 2 + 2 + 2 — оператор Лапласа. cr ¶x ¶y ¶z Это трехмерное уравнение в частных производных, реше‑ ние которого находится при определенных начальных и кра‑ евых условиях. 17 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы) Рассмотрим первую краевую задачу для одномерного урав‑ нения теплопроводности ¶U ( x , t ) ¶ 2U ( x, t ) + f ( x, t ) , 0 < x < l ¶t ¶x 2 U ( x, t ) = j ( x ) , U ( 0, t ) = m1 (t ) , U (l , t ) = m 2 (t ) . = (2.3) Эта задача является моделью распространения тепла в тон‑ ком теплоизолированном стержне длиной l с начальным распре‑ делением температуры j ( x ), и на концах которого задан режим изменения температуры, определяемый функциями m1 (t ) и m 2 (t ). Решение Решение задачи (2.3) можно получить аналитически в виде ряда Фурье и численно, используя метод конечных разно‑ стей [1]. Аналитическое решение имеет вид U ( x, t ) = m1 (t ) + 2 Ґ +е e ж np ц -з ч t и l ш sin n =1 Ґ 2 ж np ц -з ч t и l ш t x (m 2 (t ) - m1 (t )) + l ( ) npx ж 2 n ц m1 ( 0 ) - ( -1) m 2 ( 0 ) ч + jn з l и pn ш 2 ж np ц ч t l ш npx зи e sin l т0 ( ) 2 n ж ц m ў1 ( t ) - ( -1) m ў2 ( t ) ч d t, з fn ( t ) n p и ш n =1 l l pnx 2 2 pnx dx , fn (t ) = т f ( x, t ) sin dx . где jn (t ) = т j ( x ) sin l 0 l l 0 l +е e Численное решение методом конечных разностей получа‑ ют, заменяя само уравнение, а также начальные и краевые ус‑ ловия разностными аналогами. Для построения разностной схе‑ мы воспользуемся равномерной сеткой с шагом h по х и шагом τ по t: xi = h Ч i, t j = t Ч j , i = 0, 1, ј, n, j = 0, 1, , m. В этом случае l h = и tm = T . n 18 2.2. Некоторые модели математической физики В явных схемах значение сеточной функции на следующем временном слое сетки явно выражается через значения на пре‑ дыдущем слое. Это имеет место, например, при использова‑ нии шаблона (x ,t ) i (x i -1 ,t j ) j +1 (x ,t ) i j (x i +1 ,t j ) Для дифференциальных операторов в точке ( xi ,t j ) использу‑ ем разностные отношения U ( xi , t j +1 ) -U ( xi , t j ) ¶ 2U ¶U , Ю Ю t ¶t ¶x 2 U ( xi +1 , t j ) - 2 ЧU ( xi , t j ) +U ( xi -1 , t j ) Ю . h2 Разностное уравнение имеет вид или U i j +1 -U i j U i j+1 - 2 ЧU i j +U i j-1  j = + fi t h2 ж t ц U i j +1 = U i j + з 2 ч U i j+1 - 2 ЧU i j +U i j-1 + t Ч fi j , иh ш i = 1, 2, ј, n - 1, j = 0, 1, ..., m - 1; ( ) (2.4) U 0j = m1 (t j ) = m1 j , U nj = m 2 (t j ) = m 2 j , j = 0, 1, ј, m, U i0 = j ( xi ) = ji , i = 1, ј, n - 1, где fi j может определяться по формуле fi j = f ( xi , t j ). 19 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы) Рассмотрим решение первой краевой задачи с началь‑ ным условием j = x (1 - x ) и краевыми условиями вида m1 (t ) = t , m 2 (t ) = sin t . На рис. 2.5 представлен график решения, полученного из аналитического выражения с учетом в сумме первых пяти членов. 3 3 U( x , 0) 2 U( x , 1) U( x , 2) 1 U( x , 3) -1 -1 0.2 0.4 0.6 x 0.8 1 Рис. 2.5 Листинг программы численного решения. n := 9 d := 1 h := 20 d t := 0.005 m1(t ) := t m 2(t ) := sin t j := 0 f ( x,t ) := 0 n 2.2. Некоторые модели математической физики Результат решения показан на рис. 2.6. 3 3 U( 0) i 2 U( 200) i U( 400) i 1 U( 600) i -1 -1 0.2 0.4 0.6 iЧh 0.8 1 1 Рис. 2.6 21 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы) 2.3. Задачи 1. Обращением модели Мальтуса является модель Ферхюль‑ ста, в которой вводится ограничение на рост популяции, и ко‑ эффициент прироста уменьшается при приближении к мак‑ симально допустимому значению численности популяции. Модель имеет вид: ж N - N (t ) ц dN (t ) = k зз max чч N (t ) , N ( 0 ) = N 0 , N max dt и ш где N max — максимальная численность: а) исследовать стационарные состояния; б) получить ана‑ литическое и численное решение при различных значениях коэффициента k; в) провести анализ полученных результатов. 2. Получить решения в модели Лотки — Вольтерра с пере‑ менным коэффициентом a (t ) = A cos ( wt ) при различных значе‑ ниях параметров A и ω. 3. Получить решения и провести анализ модели dN (t ) = ( a1 (t ) - b12 (t ) M (t ) - b11 (t ) N (t )) N (t ) dt dM (t ) = ( a 2 (t ) - b21 (t ) N (t ) - b22 (t ) M (t )) M (t ) dt при определенных начальных условиях. 4. Получить решение в модели осаждения при увеличении массы частиц в процессе осаждения. Модель имеет вид g (rТ - r ж )V - 3cr ж U 2 (t ) dU 3yrТ dD - (rТ + xr ж )V = 0, dt D (t ) dt 4D (t ) где ψ — коэффициент неравномерности наращивания объема, а c — коэффициент сопротивления частицы в несущем слое. Рассмотреть случай D (t ) = D0 (1 + t ) . 5. Получить численное решение первой краевой задачи для ж px ц уравнения теплопроводности при j ( x ) = 0, f ( x, t ) = sin з ч . и l ш 22 3. Модели оптимизации 3.1. Условный и безусловный экстремум дифференцируемой функции В простейшей модели оптимизации используется мате‑ матический аппарат, применяемый для отыскания без‑ условного или условного экстремума функции несколь‑ ких переменных f ( x1 , x2 , ..., xn ) . Если эта функция дважды дифференцируема, то задача поиска безусловного экстремума на множестве E n сводится к решению на этом множестве си‑ ¶f = 0, i = 1, 2, ..., n (необходимые условия стемы уравнений ¶xi экстремума) и к исследованию в найденной точке (точках) на знакоопределенность матрицы из вторых производных ¶2 f , i = 1, 2, ..., n, j = 1, 2, ..., n (достаточные условия). ¶xi ¶x j Для определения точек условного экстремума при наличии ограничений в виде системы уравнений gk(x1, x2, ..., xn) = 0, ..., xn ) = 0, k = 1, 2, ..., m m < n используется либо метод исключения зави‑ симых переменных, либо метод множителей Лагранжа, в кото‑ ром задача сводится к исследованию на экстремум функции Лагранжа m Y ( x1 , x2 , ..., xn , l1 , l 2 , ..., l m ) = f ( x1 , x2 , ..., xn ) + е l k g k ( x1 , x2 , ..., xn ), k =1 23 3. Модели оптимизации т. е. решению системы ¶Y = 0, i = 1, 2, ..., n ¶xi ¶Y = g ( x1 , x2 , ..., xn ) = 0, k = 1, 2, ..., m ¶l k (необходимые условия). Достаточные условия проверяют, иссле‑ дуя знакоопределенность второго дифференциала функции Y ( x1 , x2 , ..., xn , l1 , l 2 , ..., l m ) с учетом связей на дифференциа‑ лы переменных x1 , x2 , ..., xn . Эти методы изложены во всех кур‑ сах математического анализа для высшей школы. Практическое решение системы уравнений (в общем случае трансцендентных) представляет определенные трудности при большом числе пере‑ менных и сложной структуре исследуемой функции и функций связей. Поэтому наряду с аналитическими методами поиска то‑ чек экстремума используются и численные методы, основанные на пошаговом определении последовательности точек M 1 , M 2 , ..., M n , ..., которая сходится к точке экстремума. При этом некоторые численные методы не требуют даже дифферен‑ цируемости исследуемой функции (методы поиска), некоторые требуют существования первых производных (градиентные ме‑ тоды) или высших производных (метод Ньютона — Рафсона) [5]. ПРИМЕР 5. Найдем точки экстремума функции x f ( x, y ) = -8xe - x + ye y . Аналитическим методом получаем систему x м -x -x y e xe e 8 + 8 + =0 пп x x н пe y - x e y = 0. по y Из второго уравнения следует x = y . Подставляя в первое, по‑ лучим трансцендентное уравнение 8 (1 - x ) = e1+ x . Численное ре‑ 24 3.1. Условный и безусловный экстремум дифференцируемой функции шение этого уравнения дает значение (с точностью до третье‑ го знака) x = 0.461. Листинг программы численного решения уравнения 8 (1 - x ) = e1+ x . x := 0 f ( x) := -8( 1 - x) + e ff ( x , y) := -8xЧe -x + y Чe x y 1+ x r := root ( f ( x) , x) = 0.461 ff ( 0.461 , 0.461) = -1.073. Марица из производных второго порядка x x 1 y ¶2 f ¶2 f x -x -x = + 16 e 8 xe e , = - 2 ey, 2 y ¶x ¶y ¶x y в точке ( 0.461, 0.461) имеет вид x ¶2 f x 2 y = e ¶y 2 y 3 ж 13.657 -5.895 ц ч, з и -5.895 5.695 ш т. е. точка является точкой минимума. Для численного решения этой задачи рассмотрим простей‑ ший метод покоординатного спуска. В этом методе координа‑ ты точек последовательности вычисляются следующим обра‑ зом. Пусть ( x0 , y0 ) — начальная точка и α — шаг по каждой из координат (шаг может быть разным по разным координа‑ там). Тогда координата x1 вычисляется по формуле x1 = x0 - a ¶f ( x 0 , y 0 ) , ¶x и в результате имеем точку ( x1 , y0 ) . Если выполнено условие f ( x1 , y0 ) < f ( x0 , y0 ) , то в точке ( x1 , y0 ) производится спуск по y, т. е. находим y1 = y0 - a ¶f ( x1 , y0 ) ¶y . 25 3. Модели оптимизации Если условие f ( x1 , y0 ) < f ( x0 , y0 ) не выполнено, то спуск по y проводим в точке ( x0 , y0 ) т. е. находим y1 = y0 - a ¶f ( x 0 , y 0 ) ¶y . В результате могут быть получены точки трех видов ( x1 , y0 ) , ( x0 , y1 ), ( x1, y1 ) . Полученную точку берем за исходную и про‑ должаем процедуру до тех пор, пока значение функции в най‑ денной точке будет меньше, чем в предыдущей. Если на пер‑ вом шаге не удалось уйти из исходной точки, то она приближенно уже является точкой минимума и нужно умень‑ шить размер шага для уточнения координат точки минимума. Листинг программы. x ff ( x , y) := -8xЧe x y -x ж и + y Чe dffy ( x , y) := e Чз 1 - y x y dffx( x , y) := e - 8 Чe -x Ч( 1 - x) xц ч yш fmin ( f , dfx , dfy , x0 , y0 , a , n) := x0 ¬ x0 y0 ¬ y0 U0 ¬ f ( x0 , y0) for k О 1 .. n xk ¬ xk -1 - a Чdfx( xk -1 , yk -1) F ¬ f ( xk , yk -1) UXk ¬ F if F < Uk -1 otherwise UXk ¬ Uk -1 26 xk ¬ xk -1 yk ¬ yk -1 - a Чdfy (xk , yk -1) F ¬ f ( xk , yk ) F ¬ f ( xk , yk -1) UXk ¬ F if F < Uk -1 otherwise 3.1. Условный и безусловный экстремум дифференцируемой функции UXk ¬ Uk -1 xk ¬ xk -1 yk ¬ yk -1 - a Чdfy (xk , yk -1) F ¬ f ( xk , yk ) Uk ¬ f (xk , yk ) if F < UXk otherwise yk ¬ yk -1 Uk ¬ UXk r ¬ Uk xr ¬ xk yr ¬ yk n ¬ ( r xr yr ) n q := fmin ( ff , dffx , dffy , 1 , 1 , 0.05 , 60) = ( -1.073 0.461 0.461 ) . Этот же результат можно было получить, используя проце‑ дуры пакета MathCAD. Листинг программы. x -x f ( x , y) := -8xЧe + y Чe y x := 1 y := 1 ж 0.461 ц ч и 0.461 ш q := Minimize ( f , x , y) = з f ( 0.461, 0.461) = -1.073. Заметим, что численные методы хорошо работают, если вы‑ бор начальной точки достаточно обоснован. В противном слу‑ чае результат будет неверным, что видно из представленного ниже листинга. 27 3. Модели оптимизации Листинг программы. x -x f ( x , y) := -8xЧe + y Чe y x := -1 y := -1 ж 2.824 ґ 104 ц ч. q := Minimize ( f , x , y) = з з 20 ч и -1 ґ 10 ш Ограничения на область, в которой ищется экстремум, могут носить не только характер уравнений, но и неравенств. В этом случае возникают задачи линейного, нелинейного и динами‑ ческого программирования. 3.2. Задачи линейного программирования В задачах линейного программирования ищется экстремум линейной функции (целевой функции) на множестве, ограни‑ ченном системой линейных уравнений и неравенств. Матема‑ тическая модель общей задачи линейного программирования имеет вид n f ( x1 , x2 , , xn ) = е ci xi ® min ( max ) , i =1 n еa j =1 n еa j =1 n еa j =1 28 ij ij ij x j Ј bi , i = 1, 2, , k , x j = bi , i = k , k + 1, , l , (3.1) x j і bi , i = l + 1, l + 2, , m, x j і 0, j = 1, 2, , p, p Ј n. 3.2. Задачи линейного программирования Множество точек M ( x1 , x2 , , xn ) , удовлетворяющих этим ограничениям, образуют допустимое множество. Далее будем считать, что все переменные неотрицательны, т. е. p = n. Вве‑ дением дополнительных переменных задача сводится к кано‑ ническому виду, в котором все ограничения являются равен‑ ствами [5, 6], т. е. получаем задачу n еc x i =1 n еa j =1 ij i i ® min ( max ) , x j = bi , i = 1, 2, , k , k Ј n, (3.2) x j і 0, j = 1, 2, , n. В основе методов решения этой задачи лежит следующее ут‑ верждение: если задача (3.1) имеет решение, то экстремум це‑ левой функции достигается в одной из угловых точек допусти‑ мого множества. Простейшим методом решения задачи (3.1) является графи‑ ческий метод, который с успехом может применяться в случае, когда целевая функция зависит от двух (максимум трех) пере‑ менных. Ниже будут рассмотрены примеры использования это‑ го метода при отыскании целочисленных решений. Общим методом решения задачи (3.1) является симплекс-ме‑ тод, описание которого можно найти во всех учебниках по ли‑ нейному программированию (см., например [6]). В практиче‑ ском отношении нас будут интересовать процедуры Minimize и Maximize пакета MathCAD, позволяющие решить эту задачу. Рассмотрим несколько практических задач, математической моделью которых является задача линейного программирова‑ ния в виде (3.1) или (3.2). 1. Задача о рационе Имеется n продуктов и m питательных веществ. Известны параметры: ai j — содержание (в весовых единицах) в i‑м про‑ дукте j‑го питательного вещества, b j — минимально допустимая 29 3. Модели оптимизации потребность (в единицу времени) в j‑м веществе, xi — потреб‑ ность в i‑м продукте (в единицу времени), ci — стоимость весо‑ вой единицы i‑го продукта. Требуется найти xi так, чтобы обеспечить потребность в пи‑ тательных продуктах за минимальную стоимость. Математическая модель n еc x i =1 n еa i =1 ij i i ® min, xi = b j , j = 1, 2, , m, xi і 0, i = 1, 2, , n ПРИМЕР 6. Металлургический завод выпускает три вида сплавов из двух видов сырья типа А и Б. Распределение сы‑ рья (в процентах) по трем видам сплавов представлено в табл.: Вид сырья Номер сплава А Б 1 2 3 10 % 90 % 15 % 85 % 30 % 70 % Завод имеет 6 тонн сырья типа А и 18 тонн сырья типа Б. Тон‑ на первого сплава стоит 1600 р., второго 1550 р., и третьего — 1300 р. Кроме того, известно, что спрос на сплав № 3 не пре‑ вышает 13 тонн, а сплава № 1 нужно произвести не менее пяти тонн. Сколько нужно выпустить сплавов первого, второго и тре‑ тьего типов, чтобы получить наибольшую стоимость? 30 3.2. Задачи линейного программирования Решение Математическая модель Целевая функция f (m1 , m2 , m ) = 1600m1 + 1550m2 + 1300m3 ® max . Ограничения 0.1m1 + 0.15m2 + 0.3m3 Ј 6, 0.9m1 + 0.85m2 + 0.7m3 Ј 18, m1 і 5, m2 і 0, m3 і 0, m3 Ј 13. Листинг программы. m1 := 6 m2 := 9 m3 := 7 f ( m1 , m2 , m3) := 1600 Чm1 + 1550 Чm2 + 1300 Чm3 Given m1 і 5 m2 і 0 m3 і 0 m3 Ј 13 0.1 Чm1 + 0.15 Чm2 + 0.3 Чm3 Ј 6 0.9 Чm1 + 0.85 Чm2 + 0.7 Чm3 Ј 18 ж 5 ц з ч P := Maximize ( f , m1 , m2 , m3) = 5.176 з ч и 13 ш f (P 0 , P 1 , P 2) = 3.292 ґ 10 . 4 Таким образом, максимальная стоимость произведенной продукции достигается при производстве 5 тонн первого спла‑ ва, 5.176 тонн второго и 13 тонн третьего. 2. Транспортная задача Целью транспортной задачи является составление плана пе‑ ревозок однородного груза так, чтобы стоимость перевозок была минимальной. Пусть m — число пунктов отправки груза, ai — количество единиц груза в i‑м пункте отправления (i = 1, 2, , m ) , b j — по‑ 31 3. Модели оптимизации требность в грузе в j‑м пункте приема ( j = 1, 2, , n ) , ci j — сто‑ имость перевозки из i‑го пункта отправления в j‑й пункт при‑ ема, xi j — стоимость перевозки из i -го пункта отправления в j‑й пункт приема. Математическая модель (замкнутая транспортная задача) m n е еc ij i =1 j =1 xi j ® min, m еx i =1 n еx j =1 ij ij = bj , = ai , xi j і 0, i = 1, 2, , m, j = 1, 2, , n. В замкнутой транспортной задаче суммарное количество гру‑ за в пунктах отправления равно суммарной потребности в пун‑ m n i =1 j =1 ктах приема, т. е. е ai = е b j . ПРИМЕР 7. Два цементных завода поставляют продукцию трем домостроительным комбинатам. Суточные объемы про‑ изводства цемента на заводах, потребности комбинатов и сто‑ имость перевозки цемента представлены в таблице: Заводы 1 2 32 Производство цемента (т) 40 60 Потребность в цементе Стоимость перевозки одной тонны (р.) Комб. 1 Комб. 2 Комб. 3 100 200 150 300 250 300 50 20 30 3.2. Задачи линейного программирования Листинг программы. ж 100 150 250 ц ч и 200 300 300 ш Матрица транспортных издержек p := з ж1 1 1ц ч и1 1 1ш Начальные значения x := з 1 Целевая функция f (x) := 2 е е (pi , j Чxi , j) i=0 j=0 Given ж 1 ц з xi , 0 ч з ч иi = 0 ш ж 1 ц 50 з xi , 1 ч з ч иi = 0 ш ж 2 ц з x0 , j ч з ч иj = 0 ш 40 з е е ж е 2 е з иj = 0 ц x1 , j ч ч ш ж 1 ц 20 з xi , 2 ч з ч иi = 0 ш е 30 60 x і 0 ж 20 20 0 ц 4 ч f ( P ) = 2 ґ 10 . и 30 0 30 ш P := Minimize ( f , x) = з Из полученного решения следует, что поставки с первого за‑ вода на третий комбинат и со второго завода на второй комби‑ нат отсутствуют, а минимальные издержки составляют 20000 р. m n i =1 j =1 Если е ai № е b j , то имеем открытую транспортную задачу. Рассмотрим два возможных случая (примеры 8 и 9) в задаче ПРИМЕРА 7. 33 3. Модели оптимизации ПРИМЕР 8. Избыток производства цемента. Листинг программы. ж 100 150 250 ц Матрица транспортных издержек p := з ч и 200 300 300 ш ж1 1 1ц Начальные значения x := з ч и1 1 1ш 1 Целевая функция f (x) := 2 (pi , j Чxi , j) е е i=0 j=0 Given ж 1 ц з xi , 0 ч з ч иi = 0 ш е 2 е j=0 x0 , j ж 50 з 1 е з иi = 0 ц xi , 1 ч ч ш ж 20 з 1 е з иi = 0 ц xi , 2 ч ч ш 30 ж 2 ц x1 , j ч Ј 65 x і 0 Ј 55 з ч з иj = 0 ш е ж 35 20 0 ц 4 ч f ( P ) = 1.85 ґ 10 . и 15 0 30 ш P := Minimize ( f , x) = з Стоимость перевозок уменьшилась за счет увеличения про‑ изводства цемента на первом заводе. ПРИМЕР 9. Недостаток производства цемента. Листинг программы. ж 100 150 250 ц Матрица транспортных издержек p := з ч и 200 300 300 ш ж1 1 1ц Начальные значения x := з ч и1 1 1ш 34 3.2. Задачи линейного программирования 1 Целевая функция f (x) := 2 е е (pi , j Чxi , j) i=0 j=0 Given е xi , 0 Ј 50 е xi , 1 Ј 50 2 x0 , j j=0 е xi , 2 Ј 50 i=0 i=0 i=0 е 1 1 1 ж 2 ц 40 з x1 , j ч з ч иj = 0 ш е 60 x і 0 ж 0 40 0 ц 4 ч f ( P ) = 1.9 ґ 10 . 50 10 и ш P := Minimize ( f , x) = з В данном случае весь цемент, произведенный на первом за‑ воде, поставляется на второй объект, а цемент со второго завода поставляется на первый объект (50 т.) и на третий объект (10 т.). 3. Задача о назначениях Пусть есть n типов работ и n исполнителей, ci j — стоимость работы i-го исполнителя при выполнении j‑й работы. Требу‑ ется так распределить исполнителей по работам, чтобы стои‑ мость выполнения работ была минимальной. Это чисто комбинаторная задача, в которой нужно выбрать оптимальный маршрут в таблице стоимостей: Работы Исполн. 1 2 1 2 … n c11 c21 c12 c22 … … c1n c2n n cn1 cn2 … cnn 35 3. Модели оптимизации Если выбрать по одному значению из каждой строки и каж‑ дого столбца и складывать, то получим n! чисел, из которых нуж‑ но выбрать наименьшее. При больших значениях n задача тре‑ бует работы компьютера. Для решения этой задачи можно использовать модель транс‑ портной задачи, если положить, м1, если i ‰ый исполнитель выполняет j ‰ю работу xi j = н о0, в противном слуучае. Тогда, если ci j — стоимость выполнения i‑м исполнителем j‑й работы, то математическая модель имеет вид n n е еc ij i =1 j =1 xi j ® min, n еx i =1 n еx ij j =1 n ij = 1, = 1, xi j і 0, i = 1, 2, , n, j = 1, 2, , n. n Так как е е xi j = n, то это замкнутая транспортная задача. i =1 j =1 ПРИМЕР 10. Управление строительством имеет в распо‑ ряжении 5 кранов, которые нужно распределить по пяти стро‑ ящимся объектам. Стоимость эксплуатации каждого крана на каждом объекте приведена в таблице: Кр. 1 Кр. 2 Кр. 3 Кр. 4 Кр. 5 36 Объект 1 30 20 40 90 60 Объект 2 70 40 70 70 40 Объект 3 50 40 20 30 30 Объект 4 80 50 80 80 60 Объект 5 60 70 90 100 70 3.2. Задачи линейного программирования Листинг программы. Матрица издержек эксплуатации и матрица начальных зна‑ чений ж 30 з 20 з p := з 40 з 90 з и 60 70 50 80 60 ц 40 40 50 70 20 80 70 30 80 40 30 60 ч 70 ч 90 ч x := 100 ч ч 70 ш 4 ж0 з0 з з0 з0 з и0 4 е е Целевая функция f ( x) := 0 0 0 0ц ч ч 0 0 0 0ч 0 0 0 0ч ч 0 0 0 0ш 0 0 0 0 (pi , j Чxi , j) i=0 j=0 Given 4 е 4 xi , 0 1 е i=0 i=0 4 4 е xi , 3 1 4 j=0 1 е xi , 2 1 x2 , j 1 i=0 xi , 4 1xі 0 x1 , j 1 i=0 i=0 е е 4 xi , 1 4 x0 , j 1 е 4 j=0 ж0 з0 з P := Minimize ( f , x) = з 1 з0 з и0 е 4 j=0 е 4 x3 , j j=0 1 е x4 , j 1 j=0 0 0 0 1ц ч ч 0 0 0 0 ч f ( P ) = 220. 0 1 0 0ч ч 1 0 0 0ш 0 0 1 0 37 3. Модели оптимизации Матрица Р дает распределение кранов по объектам. Из полу‑ ченных данных следует, что 1‑й кран следует отправить на 5‑й объект, 2‑й кран на 4‑й, третий на первый и т. д. Заметим, что решение не единственное. В частности, если ж0 з1 з Р = з0 з з0 з0 и 1 1 1 1ц 0 чч 0 ч, ч 0ч 0 чш то получаем то же самое значение целевой функции. ПРИМЕР 11. Пусть в таблице ПРИМЕРА 10 не будет пято‑ го объекта, т. е. число кранов превышает число объектов стро‑ ительства. Объект 1 30 20 40 90 60 Кр. 1 Кр. 2 Кр. 3 Кр. 4 Кр. 5 Объект 2 70 40 70 70 40 Объект 3 50 40 20 30 30 Объект 4 80 50 80 80 60 Листинг программы. ж 30 з 20 з p := з 40 з 90 з и 60 70 50 80 ц 40 40 70 20 70 30 40 30 4 f ( x) := 3 е е ч 50 ч 80 ч x := 80 ч ч 70 ш ж0 з0 з з0 з0 з и0 0 0 0ц ч ч 0 0 0ч 0 0 0ч ч 0 0 0ш 0 0 0 (pi , j Чxi , j) i=0 j=0 Given 38 4 е i=0 4 xi , 0 1 е i=0 4 xi , 1 1 е i=0 4 xi , 2 1 е i=0 xi , 3 1xі 0 з ч и 60 40 30 70 ш 4 f ( x) := 3 ч з и0 0 0 0ш (pi , j Чxi , j) е е 3.2. Задачи линейного программирования i=0 j=0 Given 4 е 4 xi , 0 1 i=0 i=0 3 е е 4 xi , 1 1 е i=0 3 x0 , j Ј 1 j=0 е 4 xi , 2 1 е j=0 ж1 з0 з P := Minimize ( f , x) = з 0 з0 з и0 е 1xі 0 i=0 3 x1 , j Ј 1 xi , 3 3 x2 , j Ј 1 j=0 е 3 x3 , j Ј 1 j=0 е x4 , j Ј 1 j=0 0 0 0ц ч ч 0 1 0 ч f ( P ) = 140. 0 0 0ч ч 1 0 0ш 0 0 1 Таким образом, четвертый кран вообще никуда не направ‑ ляется. В примерах 10 и 11 решения получились целочисленными в силу того, что все вершины допустимого множества имели целочис‑ ленные координаты. Однако может случиться так, что, исполь‑ зуя процедуры Miniize и Maximize, мы получим нецелочисленные решения, хотя по смыслу задачи решения должны быть целыми. ПРИМЕР 12. В доме общей площадью 5000 м 2 планируется разместить однокомнатные, двухкомнатные, трехкомнатные и четырехкомнатные квартиры площадью 50.6, 100.5, 150.8, 230.6 м 2 соответственно. По требованию заказчика одноком‑ натных и двухкомнатных квартир должно быть не менее 10 каж‑ дого типа, а трехкомнатных и четырехкомнатных — не менее 7. Прибыль застройщика с квартир разного типа составляет: Прибыль (тыс. руб.) Однокомн. Двухкомн. Трехкомн. Четырехкомн. 190.2 295.0 650.1 798.5 39 3. Модели оптимизации Требуется определить, при каком соотношении квартир раз‑ ных типов прибыль застройщика будет максимальной. Листинг программы. f ( n1 , n2 , n3 , n4) := 190.2 Чn1 + 295.0 Чn2 + 650.1 Чn3 + 798.5 Чn4 n1 := 20 n2 := 20 n3 := 20 n4 := 20 Given 50.6 Чn1 + 100.5 Чn2 + 150.8 Чn3 + 230.6 Чn4 Ј 5000 n1 і 10 n2 і 10 n3 і 7 n4 і 7 ж 10 ц з ч 10 ч 4 з P := Maximize ( f , n1 , n2 , n3 , n4) = f ( 10 , 10 , 12.432 , 7) = 1.852 ґ 10 . з 12.432 ч з ч и 7 ш Как видно из полученного результата, число трехкомнат‑ ных квартир является дробным. Возникает желание взять бли‑ жайшее целое значение и использовать его в качестве реше‑ ния. Ниже на примерах будет показано, что это может привести к неверному результату. В рамках данного пособия нет воз‑ можности рассматривать вопросы целочисленного линейного программирования, в силу чего ограничимся ссылкой, напри‑ мер [6]. После рассмотрения нескольких примеров, демонстри‑ рующих ситуации, возникающие в целочисленных задачах, мы вернемся к решению ПРИМЕРА 12. ПРИМЕР 13. Рассмотрим три двумерные задачи. 1. f ( x1 , x2 ) = -4 x1 - 3x2 ® min, 4 x1 + x2 Ј 0, 2 x1 + 3x2 Ј 8, x1 і 0, x2 і 0. Рассмотрим графическое решение этой задачи (рис. 3.1). 40 3.2. Задачи линейного программирования 4 3 2 1 1 2 3 4 Рис. 3.1 На рис. 3.1 пунктирными линиями изображены линии уров‑ ня 4 x1 - 3x2 = C . Стрелка указывает направление увеличения зна‑ чения параметра C. Одна из линий проходит через вершину до‑ пустимой области, а вторая проходит через ближайшую целую точку ( 2, 1) , которая соответствует целочисленному решению. Координаты нецелочисленной вершины находим из системы м4 x1 + x2 = 10 . В результате имеем точку ( 2.2, 1.2 ) . Тот же резуль‑ н о2 x1 + 3x2 = 8 тат можно получить, используя пакет MathCAD. Листинг программы. f1 ( x1 , x2) := -4 Чx1 - 3 Чx2 x1 := 1 x2 := 1 Given 4 Чx1 + x2 Ј 10 2 Чx1 + 3 Чx2 Ј 8 x1 і 0 x2 і 0 ж 2.2 ц ч f1 (X0 , X1) = -12.4 . и 1.2 ш X := Minimize ( f1 , x1 , x2) = з 41 3. Модели оптимизации Заметим, что значение целевой функции в точке ( 2, 1) рав‑ но –11. В данном примере целочисленное решение было дей‑ ствительно получено в ближайшей точке. 2. Рассмотрим задачу f ( x1 , x2 ) = - x1 - x2 ® min, 2 x1 + 3x2 Ј 36, 3x1 + x2 і 6, x1 і 0, x2 і 0, x1 Ј 13. Как и в предыдущем случае, минимальное значение целе‑ вой функции достигается в точке с нецелочисленными коор‑ м2 x + 3x2 = 36 10 ц ж , которые находим из системы н 1 динатами з13, ч 3 ш и о x1 = 13 или используя программу. Листинг программы. f2 ( x1 , x2) := -x1 - x2 x1 := 1 x2 := 1 Given 2 Чx1 + 3x2 Ј 36 3 Чx1 + x2 і 6 x1 і 0 x2 і 0 x1 Ј 13 ж 13 ц ч f2 (X0 , X1) = -16.333 и 3.333 ш X := Minimize ( f2 , x1 , x2) = з На рис. 3.2 представлены допустимая область — а) и реше‑ ние — б). а) б) 12 6 5 4 3 6 12 Рис. 3.2 42 9 10 11 12 13 3.2. Задачи линейного программирования Видно, что в данном случае решение не единственное. В точ‑ ках (13, 3) и (12, 4 ) значение целевой функции равно –16. 3. Рассмотрим задачу f ( x1 , x2 ) = x1 - 20 x2 ® min, - x1 + 5x2 Ј 19, 4 x1 + 2 x2 Ј 29, x1 і 0, x2 і 0. 17 ц ж 10 Точка минимума имеет координаты з 4 , 4 ч . 22 ш и 22 Листинг программы. f3 ( x1 , x2) := x1 - 20x2 x1 := 1 x2 := 1 Given -x1 + 5x2 Ј 19 4 Чx1 + 2x2 Ј 29 x1 і 0 x2 і 0 ж 4.864 ц ч и 4.773 ш X := Minimize ( f3 , x1 , x2) = з f3 (X0 , X1) = -90.591 . Графически имеем ситуацию, показанную на рис. 3.3. 4 3 2 1 1 2 3 4 5 6 Рис. 3.3 43 3. Модели оптимизации Таким образом, целочисленное решение получаем в точке (1, 4 ), которая не принадлежит ближайшей окрестности точки вершины. Значение целевой функции в этой точке равно –70. Заметим, что задачи целочисленного программирования, в принципе, можно решать простым перебором целых точек в допустимой области. В частности в задаче 2 листинг програм‑ мы перебора имеет вид x1 := 8 .. 13 x2 := 0 .. 7 Fx1 -8 , x2 := -x1 - x2 if 2x1 + 3x2 Ј 36 Щ ( 3x1 + x2) і 6 Щ x1 Ј 13 0 otherwise ж -8 з з -9 з -10 F=з з -11 з -12 з -13 и -9 -10 -11 -12 -13 -14 0 ц ч -10 -11 -12 -13 -14 -15 0 ч -11 -12 -13 -14 -15 0 0 ч -12 -13 -14 -15 -13 -14 -15 -16 -14 -15 -16 ч 0ч min ( F) = -16 . 0ч ч 0ш Из полученного результата следует, что минимальное зна‑ чение –16 достигается в точках (12, 4 ) и (13, 3) . Вернемся к рассмотрению задачи ПРИМЕРА 12. Целочис‑ ленное решение найдем перебором. Листинг программы. f ( x1 , x2 , x3 , x4) := 190.2 Чx1 + 295.0 Чx2 + 650.1 Чx3 + 798.5 Чx4 FF := f0 ¬ f ( 10 , 10 , 7 , 7) for n1 О 10 .. 98 for n2 О 10 .. 49 for n3 О 7 .. 33 for n4 О 7 .. 21 if 50.6 Чn1 + 100.5 Чn2 + 150.8 Чn3 + 230.6 Чn4 Ј 5000 44 F ¬ f ( n1 , n2 , n3 , n4) if F і f0 for n1 О 10 .. 98 for n2 О 10 .. 49 for n3 О 7 .. 33 3.3. Нелинейное программирование for n4 О 7 .. 21 if 50.6 Чn1 + 100.5 Чn2 + 150.8 Чn3 + 230.6 Чn4 Ј 5000 F ¬ f ( n1 , n2 , n3 , n4) if F і f0 f0 ¬ F p ¬ ( F n1 n2 n3 n4 ) ( p FF = 1.843 ґ 10 4 11 10 12 7 ) Максимальное значение целевой функции равно 18430 р. Точка, в которой достигается это значение, имеет координаты (11, 10, 12, 7) . Заметим, что вопрос об единственности решения нужно исследовать дополнительно. Изменяя начальное значе‑ ние параметра n2 с 10 на 11, убеждаемся, что полученная точка единственная. 3.3. Нелинейное программирование В пособии будут указаны только некоторые задачи нелиней‑ ного программирования, которые можно решить, используя возможности пакета МathCAD. Это задачи дробно-линейного программирования, квадратичного программирования и неко‑ торые задачи выпуклого программирования. 1. Дробно-линейное программирование Типичными задачами дробно-линейного программирова‑ ния являются задачи на минимизацию себестоимости про‑ дукции и максимизацию производительности. Математиче‑ ская модель задачи дробно-линейного программирования имеет вид 45 3. Модели оптимизации n f ( x1 , x2 , , xn ) = еa x i + a9 еb x + b9 i i =1 n i i =1 n еc ij j =1 n еc j =1 ij i ® min, (3.3) x j = di , i = 1, 2, , k , x j Ј di , i = k + 1, 2, , m, xi і 0, i = 1, 2, , n. Считается, что знаменатель целевой функции положителен на допустимом множестве. xi 1 , i = 1, 2, , n, x0* = n Заменой xi* = n задача + b x b b x + b еi i 0 еi i 0 i =1 i =1 (3.3) сводится к задаче линейного программирования n f ( x1 , x2 , , xn ) = е ai xi* ® min, i =1 n еc j =1 n еc j =1 ij ij * j (3.4) * x - di x = 0, i = 1, 2, , k , x *j - di x0* Ј 0, i = k + 1, 2, , m, xi* і 0, i = 1, 2, , n. ПРИМЕР 14. Три портовых крана разгружают сухогруз, на котором находятся четыре типа различных грузов в коли‑ честве 250, 310, 400 и 600 тонн. Прибыль, получаемая портом (руб.) за разгрузку тонны разных грузов разными кранами (ai j ) приведена в таблице: Кран 1 Кран 2 Кран 3 46 Груз 1 300 250 250 Груз 2 250 250 200 Груз 3 400 300 300 Груз 4 350 350 250 3.3. Нелинейное программирование Затраты времени (ч.) на разгрузку тонны груза (bi j ) приведе‑ ны ниже: Груз 1 0.018 0.030 0.020 Кран 1 Кран 2 Кран 3 Груз 2 0.015 0.020 0.020 Груз 3 0.020 0.025 0.024 Груз 4 0.022 0.030 0.025 Кроме того, известно, что первый кран за все время рабо‑ ты может выгрузить не более 450 тонн груза, второй — не бо‑ лее 550. Составить план работ так, чтобы рентабельность раз‑ грузки одной тонны груза была максимальной. Математическая модель Пусть xi j — количество груза j‑го типа, разгруженного i‑м краном. Целевая функция имеет вид 3 4 ее ai j xi j ее bi j xi j i =1 j =1 3 4 i =1 j =1 ® max, где в числителе стоит суммарная прибыль от работы всех кра‑ нов, а в знаменателе суммарное время работы. Ограничения имеют вид 3 е хi 1 = 250, i =1 4 ех j =1 1j 3 е хi 2 = 310, Ј 450, i =1 4 ех j =1 2j 3 е хi 3 = 400, i =1 Ј 550, 4 ех j =1 3j 3 ех i =1 Ј 600, i4 = 600, хi j і 0. 47 3. Модели оптимизации Листинг программы. a := ж 300 250 400 350 ц з 250 250 300 350 ч b := з ч и 250 150 300 250 ш ж 0.018 0.015 0.020 0.022 ц з 0.030 0.020 0.025 0.030 ч ч з и 0.020 0.020 0.024 0.025 ш ж 100 100 150 100 ц x := з 150 120 180 120 ч f ( x) := з ч и 125 190 100 170 ш 3 й 2 щ к (xi , j Чai , j)ъ к ъ лi = 0 j = 0 ы е е 2 3 е е (xi , j Чbi , j) ж ц i=0 j=0 Given ж 2 ц з xi , 0 ч з ч иi = 0 ш е 3 е j=0 ж 250 з x0 , j Ј 450 2 е з иi = 0 ц xi , 1 ч ч ш j=0 2 е з иi = 0 xi , 2 ч ч ш ж 400 з 2 е з иi = 0 ц xi , 3 ч ч ш 600 3 3 е 310 з x1 , j Ј 550 е x2 , j Ј 600 x і 0 j=0 ж 0 0 400.001 50 ц з 0 310 ч f ( r ) = 1.359 ґ 104. r := Maximize ( f , x) = 240 з ч 310 ш и 250 0 Таким образом, задача дробно-линейного программирова‑ ния может быть решена с использованием стандартных проце‑ дур пакета MathCAD. 2. Квадратичное программирование Пусть целевая функция квадратичная, а ограничения име‑ ют линейный характер. Типичной задачей является задача о ре‑ жиме работы нескольких производителей однородной продук‑ ции, снабжающих одного потребителя. 48 3.3. Нелинейное программирование ЗАДАЧА. Пусть есть n производителей и xi — количество продукции, которое способен произвести i -й производитель, причем ai Ј xi Ј bi . P — потребность в продукции в пункте по‑ n требления, т. е. е xi = P . Затраты на производство xi продукции i =1 есть квадратичная на [a, b] функция f ( xi ) = ai xi2 + bi xi + ci . Требу‑ ется минимизировать суммарные затраты на производство дан‑ ного объема потребления. Математическая модель n е f ( x ) ® min, i i =1 n еx i =1 i = P , ai Ј xi Ј bi . ПРИМЕР 15. Четыре электростанции снабжают электриче‑ ством один город. Потребность города в электричестве состав‑ ляет 70 МВт. Мощности электростанций представлены в та‑ блице: Электростанции Мощности Затраты (на x продукции) 1 13–15 2 16–20 3 14–20 4 18–24 x 2 - 24 x + 160 x 2 - 30 x + 240 x 2 - 28x + 300 x 2 - 36 x + 334 Решение Листинг программы. 2 2 2 f1 ( x) := x - 24 Чx + 160 f2 ( x) := x - 30 Чx + 240 f3 ( x) := x - 28 Чx + 300 2 f4 ( x) := x - 36 Чx + 334 f ( x1 , x2 , x3 , x4) := f1 ( x1) + f2 ( x2) + f3 ( x3) + f4 ( x4) x1 := 13 x2 := 15 x3 := 14 x4 := 14 Given x1 і 13 x2 і 16 x3 і 14 x4 і 18 x1 Ј 15 x2 Ј 20 x3 Ј 20 x4 Ј 24 ж 14.75 ц з ч 17.75 49 2 2 2 f1 ( x) := x - 24 Чx + 160 f2 ( x) := x - 30 Чx + 240 f3 ( x) := x - 28 Чx + 300 2 f4 ( x) оптимизации := x - 36 Чx + 334 f ( x1 , x2 , x3 , x4) := f1 ( x1) + f2 ( x2) + f3 ( x3) + f4 ( x4) 3. Модели x1 := 13 x2 := 15 x3 := 14 x4 := 14 Given x1 і 13 x2 і 16 x3 і 14 x4 і 18 x1 Ј 15 x2 Ј 20 x3 Ј 20 x4 Ј 24 x1 + x2 + x3 + x4 ж 14.75 ц з ч 17.75 ч з 70 r := Minimize ( f , x1 , x2 , x3 , x4) = з 16.75 ч з ч и 20.75 ш f ( 14.75, 17.75, 16.75, 20.75) = 175.25 . Таким образом, распределение мощностей имеет вид: 14.75, 17.75, 16.75, 20.75 для 1‑й, 2‑й, 3‑й и 4‑й электростанции соот‑ ветственно. Общая задача квадратичного программирования имеет вид n n ее i =1 j =1 n ai j xi x j + е bi xi + b0 ® min, i =1 n еc x i =1 i i Ј di , i = 1, 2, , n, xi і 0, i = 1, 2, , n, n n где е е ai j xi x j — положительно или отрицательно полуопре‑ i =1 j =1 деленная квадратичная форма. Если число переменных не превышает трех, то для задачи квадратичного программирования можно использовать графи‑ ческий метод. Продемонстрируем это на примере. ПРИМЕР 16. Решить задачу f ( x1 , x2 ) = x12 - 2 x1 - x2 ® min, 2 x1 + 3x2 Ј 6, 2 x1 + x2 Ј 4, x1 і 0, x2 і 0. Допустимая область и линии уровня целевой функции (па‑ 2 раболы ( x1 - 1) - x2 - 1 = C ) представлены на рис. 3.4. 50 3.3. Нелинейное программирование 3 2 1 1 2 3 Рис. 3.4 Координаты точки касания и значения параметра C на‑ 2 пм x - 2 x1 - x2 = C . Выражая из второго урав‑ ходим из системы н 1 по2 x1 + 3x2 = 6 нения x2 и подставляя в первое, получаем уравнение 4 x12 - x1 - 2 - C = 0. Из условия равенства нулю дискриминанта 3 4 2 x12 - 2 x1 - x2 = C имеем C = -2 и x1 = . Таким образом, коорди‑ 9 3 ж 2 14 ц наты точки касания з , ч. и3 9 ш Листинг программы. 2 f ( x1 , x2) := x1 - 2 Чx1 - x2 Given 2x1 + 3x2 Ј 6 2x1 + x2 Ј 4 x1 і 0 x2 і 0 ж 0.667 ц ч. и 1.556 ш r := Minimize ( f , x1 , x2) = з 51 3. Модели оптимизации 3.4. Выпуклое программирование Задачей выпуклого программирования является минимиза‑ ция выпуклой функции в случае, когда допустимое множество является выпуклым. Рассмотренные выше задачи линейного и квадратичного программирования являются частными слу‑ чаями выпуклого программирования. Теория задач выпуклого программирования излагается в [5]. Здесь мы рассмотрим толь‑ ко один пример использования процедур пакета MathCAD для решения подобной задачи. ПРИМЕР 17. Рассмотрим задачу f ( x1 , x2 , х3 ) = x12 - 2 x1 + х22 - 4 x2 - х3 + 5 ® min, x12 + x22 + х32 Ј 4, x1 і 0, x2 і 0, x3 і 0. Геометрическое представление задачи дает, очевидно, един‑ ственное решение, так как поверхности уровня есть параболо‑ иды с вершиной, координаты которой x0 = 1, y0 = 2, а допусти‑ мая область есть часть шара с радиусом 2. Листинг программы. 2 2 f ( x1 , x2 , x3) := x1 - 2x1 + x2 - 4 Чx2 - x3 + 5 x1 := 1 x2 := 1 x3 := 1 Given 2 2 2 x1 + x2 + x3 Ј 4 x1 і 0 x2 і 0 x3 і 0 ж 0.709 ц r := Minimize ( f , x1 , x2 , x3) = з 1.417 ч f ( r 0 , r 1 , r 2) = -0.796 . з ч и 1.221 ш Рассмотрим случай, когда целевая функция не является вы‑ пуклой в допустимой области. 52 3.4. Задачи ПРИМЕР 18. Найти решение задачи ( ) f ( x1 , x2 , х3 ) = cos 2 x12 + 2 х22 - х3 + 1 ® min, x + x + х Ј 4, x1 і 0, x2 і 0, x3 і 0. 2 1 2 2 2 3 Листинг программы. 2 2 f ( x1 , x2 , x3) := cos йл( 2x1) + 2x2 щы - x3 + 1 x1 := 1 x2 := 1 x3 := 1 Given 2 2 2 x1 + x2 + x3 Ј 4 x1 і 0 x2 і 0 x3 і 0 ж -7.312 ґ 104 ц з ч з 5 r := Minimize ( f , x1 , x2 , x3) = -6.994 ґ 10 ч . з ч з 1.282 ґ 107 ч и ш В результате получены недостоверные данные. Таким обра‑ зом, некритическое применение стандартных процедур может привести к результату, который не соответствует действитель‑ ной ситуации, или вообще давать неверный результат. 3.4. Задачи 1. Найти экстремумы функции f ( x, y ) = ln ( x 2 + y 2 + 1) при ус‑ ловии e x + y - xy = 0. 2. Рукопись, которая содержит 386,5 страниц, распределена между тремя наборщиками. Время, за которое первый набор‑ щик набирает одну страницу, равно 12 минут, второй — 15 мин, третий — 18 мин. Рукопись нужно набрать не более чем за пять смен (40 часов). Кроме того, известно, что каждый наборщик работает в смену не менее пяти часов и не более семи, и стои‑ 53 3. Модели оптимизации мость часа работы первого наборщика равна 359 рублей, второ‑ го — 400 рублей и третьего — 375 рублей. Как следует распре‑ делить рукопись между наборщиками, чтобы стоимость набора была наименьшей? 3. Найти целочисленное решение задачи 2 при условии, что рукопись имеет 380 страниц. 4. В трех складах находится суточный запас деталей типа А и В в количестве, указанном в таблице: 1‑й склад 20 30 Детали А Детали В 2‑й склад 10 20 3‑й склад 12 6 Детали со склада поставляются в два цеха комплектами. С первого завода комплект содержит две детали А и три дета‑ ли В. С второго склада — одну деталь А и две детали В. С третье‑ го склада — две детали А и одну деталь В. Затраты на поставку одного комплекта из разных складов в разные цехи приведе‑ ны в таблице: 1‑й склад 50 30 Цех 1 Цех 2 2‑й склад 30 10 3‑й склад 20 15 Суточная потребность в деталях в цехах составляет: Детали А 32 10 Цех 1 Цех 2 Детали В 40 16 Найти количество комплектов, поставляемых из разных складов в разные цеха при минимальных транспортных рас‑ ходах. 5. Найти решение задачи f ( x1 , x2 , х3 ) = x12 + х22 - х32 + 1 ® min, ( x1 - 2) + ( x2 - 1) 2 54 2 + х32 Ј 4, x1 і 0, x2 і 0, x3 і 0. 4. Дискретно-детерминированные схемы (F‑схемы) 4.1. Конечные автоматы Д искретно-детерминированные схемы описывают ра‑ боту конечных автоматов. Особенностью работы ко‑ нечных автоматов является дискретное время. Чтобы описать математическую модель конечного автомата, необ‑ ходимо задать следующие множества: X — множество вход‑ ных сигналов; Y — множество выходных сигналов; Z — мно‑ жество внутренних состояний. Кроме этого необходимо задать еще две функции: функцию переходов и функцию выходов. Функция переходов j ( x, z ) определяет внутреннее состоя‑ ние z О Z , в которое переходит система, находящаяся в состоя‑ нии z под воздействием внешнего сигнала x. Функция выходов y ( x, z ) определяет выходной сигнал y ОY , который генерирует система, находящаяся в состоянии z под воздействием внеш‑ него сигнала x. Математическая модель есть совокупность F = X , Y , Z , j ( x, z ) , y ( x, z ) . Процесс работы автомата происходит в дискретные момен‑ ты времени t : t0 , t1 ,  . Моменты времени удобно обозначать цифрами, а именно t : 0, 1, . Считается, что в начальный мо‑ мент времени t0 = 0 состояние системы было z (t0 ) = z ( 0 ) = z0 О Z . Таким образом, процесс работы имеет вид: 55 4. Дискретно-детерминированные схемы (F‑схемы) 1) t = 0, z ( 0 ) = z0 О Z , 2) в момент времени t = i, i і 0 поступает входной сигнал x (i ) О X , 3) генерируется выходной сигнал y (i ) = y ( x (i ) , z (i )) , i = 0, 1, 2, ј , 4) автомат переходит в состояние z (i + 1) = j ( x (i ) , z (i )) , i = 0, 1, 2, ј . Автомат, работа которого описывается набором равенств 3 и 4, называется автоматом первого рода (автоматом Мили). Возможны и другие способы задания работы автомата. Напри‑ мер, если y (i ) = y ( x (i - 1) , z (i )) , i = 1, 2, ј , то это автомат с за‑ паздыванием, автомат второго рода. Возможны частные случаи. Например, если функция выходов зависит только от внутрен‑ них состояний y (i ) = y ( z (i )) , i = 1, 2, ј . Такой автомат назы‑ вается автоматом Мура. Функции перехода и выхода можно задать аналитически, та‑ блично и с помощью графов. Рассмотрим табличное задание. Пусть X ( x1 , x2 , , xk ) , Y ( y1 , y2 , , ym ) , Z ( z1 , z 2 , , zl ) . Таблица переходов имеет вид: xi zj z1 z2 … zl x1 x2 φ(x1, z1) φ(x2, z1) φ(x1, z2) φ(x2, z2) … … φ(x1, zl) φ(x2, zl) xk φ(xk, z1) φ(xk, z2) … φ(xk, zl) Таблица выходов xi 56 zj z1 z2 … zl x1 x2 ψ(x1, z1) ψ(x2, z1) ψ(x1, z2) ψ(x2, z2) … … ψ(x1, zl) ψ(x2, zl) xk ψ(xk, z1) ψ(xk, z2) … ψ(xk, zl) 4.1. Конечные автоматы ПРИМЕР 19. Распределительное устройство работает в трех режимах: 1) спящий режим (0); 2) режим неполной актива‑ ции (1); 3) режим полной активации (2). Если устройство на‑ ходилось в режиме 1 или 2 и получило сигнал, частота которого менее 100 кГц, то оно переходит в спящий режим (или остается в нем). Если же устройство находилось в режиме 2, то при по‑ ступлении сигнала с частотой менее 100 кГц оно переходит в ре‑ жим неполной активации. Если частота поступившего сигна‑ ла находится в пределах 100–1000 кГц, то устройство переходит в режим неполной активации (или остается в нем), независимо от того, в каком оно было режиме. При поступлении сигнала с частотой свыше 1000 кГц устройство переходит в режим пол‑ ной активации (или остается в нем) из второго и третьего режи‑ мов, а из спящего режима переходит в режим неполной акти‑ вации. Устройство, перешедшее в режим неполной активации из спящего режима, передает на пульт диспетчера сигнал пер‑ вого типа (1), оставшееся в этом режиме устройство передает сигнал второго типа (2), а перешедшее из режима полной ак‑ тивации — сигнал третьего типа (3). Если устройство перешло в режим полной активации из режима 1, то на пульт передает‑ ся сигнал типа 4, а если осталось в режиме 2, то сигнал типа 5. При переходе в спящий режим сигнал на пульт не передается, т. е. передается сигнал типа 0. Описать множества входных и выходных сигналов, множе‑ ство внутренних состояний, таблицы переходов и выходов и по‑ лучить протокол работы при некоторой последовательности входных сигналов и данном начальном состоянии. Решение Множество внутренних состояний: Z0 — спящий режим (0); Z1 — режим неполной активации (1); Z2 — режим полной ак‑ тивации (2). 57 4. Дискретно-детерминированные схемы (F‑схемы) Множество входных сигналов: X0 — частота менее 100 кГц (0); X1 — частота в пределах 100–1000 кГц (1); X3 — частота свы‑ ше 1000 кГц (2). Множество выходных сигналов: Y0 — нет сигнала (0); Y1 — сигнал первого типа (1); Y2 — сигнал второго типа (2), Y3 — сиг‑ нал третьего типа (3), Y4 — сигнал четвертого типа (4), Y5 — сиг‑ нал второго типа (5). Таблица переходов X0 X1 X2 Z0 Z0 Z1 Z1 Z1 Z0 Z1 Z2 Z2 Z1 Z1 Z2 Таблицу выходов в этом случае определить нельзя. Функ‑ цию выходов определим аналитически: м1, z (t ) = 0 Щ z (t + 1) = 1 п п2, z (t ) = 1 Щ z (t + 1) = 1 пп3, z (t ) = 2 Щ z (t + 1) = 1 y (t ) = н . п4, z (t ) = 1 Щ z (t + 1) = 2 п5, z (t ) = 2 Щ z (t + 1) = 2 п по0, z (t + 1) = 0 Для описания процесса работы автомата используем задание функций перехода и выхода в программе MаthCAD. Листинг программы. Внутренние состояния z0 := 0 z1 := 1 z2 := 2 Матрица переходов 58 j ж z0 z0 z1 ц з ч := з z1 z1 z1 ч зz z z ч и 1 2 2ш j ж0 0 1ц = з1 1 1ч з ч и1 2 2ш 4.1. Конечные автоматы Функция переходов Zin ( n , Z0 , Xin) := j ¬ Z0 i ¬ Xin0 for t О 0 .. n - 1 zzt + 1 ¬ j i , j i ¬ Xint + 1 j ¬ zzt + 1 zz0 ¬ Z0 zz Входной поток Xin := ( 1 0 1 2 1 0 ) Функция выходов Yout ( n , Z0 , Xin) := T j ¬ Z0 i ¬ Xin0 for t О 0 .. n - 1 zzt + 1 ¬ j i , j if zzt + 1 1 yyt ¬ 1 if zzt yyt ¬ 2 if zzt 1 yyt ¬ 3 if zzt 2 if zzt + 1 2 yyt ¬ 4 if zzt 1 yyt ¬ 5 if zzt 2 i ¬ Xint + 1 j ¬ zzt + 1 yyt ¬ 0 if zzt + 1 yy 59 4. Дискретно-детерминированные схемы (F‑схемы) ж1ц з ч з1ч з0ч Zin ( 5 , z1 , Xin) = з ч з1ч з2ч з1ч и ш ж1ц з0ч з ч Yout (5 , z1 , Xin) = з 1 ч . з4ч з ч и3ш Таким образом, работа автомата при заданном входном по‑ токе и начальном состоянии z1 = 1 имеет вид: Номер такта X Z Y 1 1 1 1 1 2 1 1 3 2 1 4 4 1 2 3 5 1 4.2. Клеточные автоматы Клеточным автоматом является сеть из конечных автома‑ тов, внутренние состояния которых изменяются в дискретные моменты времени в зависимости от того, в каком исходном со‑ стоянии находится каждый элемент сети и в каком состоянии находятся элементы его ближайшего (возможно, не только бли‑ жайшего) окружения. Если сеть состоит из N элементов и z j (tn ) — состояние j‑го элемента в момент времени tn , то закон перехода может иметь ж ц вид z j (tn +1 ) = F з z j (tn ) , е zk (tn ) ч , т. е. состояние j‑го элемента k Оокр. j и ш в момент времени tn+1 является функцией его состояния в мо‑ мент tn и суммы состояний элементов его окружения. 60 4.2. Клеточные автоматы Закон перехода может быть задан и иначе. Например, если ввести еще и переменную z j (tn-1 ) , то состояние клеточного ав‑ томата будет обладать памятью. Кроме того, число и вид учиты‑ ваемых ближайших соседей может быть различным для разных элементов (нерегулярная сеть). В регулярной сети элементы за‑ нимают узлы правильной решетки. В частности, в квадратной решетке ближайших соседей можно выбрать, как это указано на рис. 4.1, двумя способами. а) б)     Рис. 4.1 Если каждый элемент клеточного автомата имеет q соседей и может находиться в r состояниях, то общее число возможных правил перехода, т. е. различных клеточных автоматов, равно 2 r qr . Действительно, если пронумеровать все внутренние состо‑ яния элементов z1 , z2 , , zr и присвоить им значения 1, 2, , r , то q Ч r есть число возможных различных значений суммы е k Оокр. j zk (tn ) и q Ч r 2 — число возможных комбинаций аргументов ж ц функции F з z j (tn ) , е zk (tn ) ч. Так как величина z j (tn+1 ) может k Оокр. j и ш принимать r значений, то общее число возможных способов 2 перехода равно r qr . Для одномерной сети (линейная сеть) с дву‑ мя соседями и двумя состояниями каждого элемента r = 2, q = 2 общее число клеточных автоматов равно 28 = 256. По динамике работы клеточные автоматы разделяют на че‑ тыре типа. 1. Клеточный автомат первого типа после конечного числа шагов переходит в однородное состояние, в котором все z j (tn ) одинаковы и не зависят от времени. 61 4. Дискретно-детерминированные схемы (F‑схемы) ПРИМЕР 20. Рассмотрим цепь из 16‑ти элементов, каждый из которых имеет два состояния z0 = 0, z1 = 1. Правило перехода имеет вид м1, z (t ) = z j (tn ) = z j +1 (tn ) = 1 . z j (tn +1 ) = н j -1 n о0, в противном случае Пусть задано начальное состояние , где серые квадраты соответствуют состоянию z1 = 1. Для крайних элементов цепи введем замыкание области, т. е. будем считать, что для крайнего левого элемента соседом сле‑ ва будет крайний правый элемент, и наоборот. Тогда на пер‑ вом такте получаем . На втором такте имеем однородное состояние. 2. Клеточные автоматы второго типа создают локализован‑ ные простые структуры. ПРИМЕР 21. Пусть цепь состоит из 16‑ти элементов и пра‑ вило переходов имеет вид м1, z j (tn ) = 1, z j -1 (tn ) + z j +1 (tn ) = 0 или z j -1 (tn ) + z j +1 (tn ) = 1 п z j (tn +1 ) = н1, z j (tn ) = 0, z j -1 (tn ) + z j +1 (tn ) = 2 п0, в противном случае. о Начальное состояние такое же, как в ПРИМЕРЕ 20. До ше‑ стого такта получаем 62 4.2. Клеточные автоматы 1‑й такт 2‑й такт 3‑й такт 4‑й такт 5‑й такт 6‑й такт , , , , , . Получаем периодическое повторение структур, возникаю‑ щих на 4‑м и 5‑м тактах. 3. В клеточных автоматах третьего и четвертого типов воз‑ можны ситуации, когда в процессе работы теряется зависи‑ мость от начальных условий или, наоборот, последовательно‑ сти состояний при различных начальных условиях различны. ПРИМЕР 22. Элементы клеточного автомата расположены в узлах квадратной решетки 8×8, и каждый элемент имеет два состояния: 0 и 1. Правило переходов имеет вид k =1 l =1 м п1, zi j (tn ) = 0, е е zi + k j +l (tn ) = 2 п k =-1 l =-1 z j (tn +1 ) = н , k =1 l =1 п0, z (t ) = 1, е е zi +k j +l (tn ) - zi j (tn ) і 3 по i j n k =-1 l =-1 63 4. Дискретно-детерминированные схемы (F‑схемы) т. е. если элемент был в состоянии 1, то в состояние 0 он перехо‑ дит, если в его ближайшем окружении, которое имеет 8 элемен‑ тов (ситуация б) на рис. 4.1), есть два элемента в состоянии 1. Если элемент был в состоянии 0, то в состояние 1 он перехо‑ дит при наличии в его окружении трех и более элементов в со‑ стоянии 1. В остальных случаях он не изменяет свое состояние. Начальное состояние — четыре активных элемента в центре квадрата. На рис. 4.2 показаны ситуации, возникающие на 1‑м, 2‑м и 3‑м тактах. Начальное состояние 1‑й такт     2-й такт 3-й такт     Рис. 4.2 Клеточные автоматы используются при моделировании ак‑ тивных сред. В частности, при моделировании гидродинамиче‑ ских течений, в задачах обработки информации, при моделиро‑ 64 4.3. Задачи вании химических реакций в растворах, процессов деформации и разрушения в неоднородных средах и некоторых других за‑ дачах. 4.3. Задачи 1. Автомат продает фруктовый сок в упаковках по 200 и 400 мл., по цене 200 и 360 рублей соответственно. Автомат принимает только купюры в 50 и 100 рублей. Описать множества входных и выходных сигналов, внутренних состояний, а также функции переходов и выходов и получить протокол работы при некоторой последовательности входных сигналов и данном на‑ чальном состоянии. 2. На перекрестке двух улиц с односторонним движением находится светофор, который может находиться только в двух состояниях: z0 — зеленый для первой улицы (красный для второй), z1 — красный для первой улицы (зеленый для второй). Каждый такт на перекресток поступают машины по каждой дороге и проезжают (при зеленом свете) по n машин (или все ма‑ шины, если их число меньше или равно n). Светофор переключа‑ ется на зеленый для улицы, на которой на данном такте больше машин. Получить протокол работы при некоторой последова‑ тельности вводных сигналов и данном начальном состоянии. 3. Рассмотреть задачу № 2 в случае, когда светофор имеет обычный вид, т. е. имеет и желтый сигнал, который подается только в течение одного такта. 4. В ПРИМЕРЕ 22 использовать циклические граничные ус‑ ловия и получить картину на 6‑м такте. 5. Описать динамику клеточного автомата 8×8 до 6‑го так‑ та, элементы которого могут находиться в трех состояниях: 0, 1, 2 — и имеют четыре соседа в ближайшем окружении. Пра‑ вило переходов имеет вид 65 4. Дискретно-детерминированные схемы (F‑схемы) м2, zi j (tn ) = 0, в окрестности есть хотя бы один элемен нт в состоянии 2 п z j (tn +1 ) = н2, zi j (tn ) = 1, в окрестности не менее трех элементов в состоянии 2 п о0, zi j (tn ) = 2 в окрестности ровно два элемента в состоянии 1. Обозначим zi j = 0 ; zi j = 1 ; zi j = 2 Рассмотреть начальные состояния а) б)     66 . 5. Дисктетно-стохастические схемы (Р‑схемы) 5.1. Вероятностные автоматы P‑схемы описывают работу вероятностных автоматов. Математическая модель Напомним, что конечный автомат детерминированного типа есть совокупность F = Z , X , Y , j, y , где Z = {z1 , z2 , , zl } — множество внутренних состояний, X = {x1 , x2 , , xk } — мно‑ жество входных сигналов, Y = {y1 , y2 , , ym } — множество выходных сигналов, j, y — функция переходов и функция выходов соответственно. Для описания вероятностного ав‑ томата рассмотрим множество пар вида XP = {( xi , z j )}, где i = 1, 2, , k , j = 1, 2, , l и множество пар вида YP = {(zi , y j )}, где i = 1, 2, , l , j = 1, 2, , m. Каждому элементу множе‑ ства XP поставим в соответствие некоторый закон распределе‑ ния на множестве пар вида YP так, что (z1, y1) (xq, zr) l ие i =1 p11 m еp j =1 ij (z1, y2) p12 … … (zl, ym–1) plm–1 (zl, ym) plm = 1. pi j трактуется как вероятность перехода автома‑ та в состояние zi с генерацией выходного сигнала y j при усло‑ вии, что автомат, находившийся в состоянии zr , получил внеш‑ 67 5. Дисктетно-стохастические схемы (Р‑схемы) ний сигнал xq . Число таких законов равно k ґ l — числу элементов в множестве XP. Множество этих законов обозна‑ чим YP. Каждый из заданных законов распределения случай‑ ного вектора ( Z ,Y ) задает закон распределения на множествах m l j =1 i =1 Z и Y по формулам p ( zi ) = е pi j и p ( y j ) = е pi j . Если pi j = p ( zi ) p ( y j ) , что соответствует статистической неза‑ висимости случайных величин Z и Y, то вероятностный авто‑ мат называется автоматом Мили. Возможна ситуация, когда выходной сигнал зависит только от внутреннего состояния, в котором находится автомат (при любом входном сигнале). В этом случае закон распределения на множестве YP ставится в соответствие не элементу ( xq , zr ) , а элементу zr (сравните с конечным детерминированным авто‑ матом Мура). Как и в конечных детерминированных автоматах, время в ве‑ роятностных автоматах дискретно и работа происходит по так‑ там. Моделью вероятностного автомата является совокупность P = Z , X , Y , YP . Для математического описания процесса работы вероятност‑ ного автомата используются цепи Маркова. Цепи Маркова. Пусть система может находиться в k различ‑ ных состояниях. Рассмотрим последовательность испытаний в моменты t1 , t2 , , tn ,  , в которой вероятность перехода в j‑е состояние на s‑м шаге зависит только от того, в каком система находилась на s-1‑м шаге. Если вероятность переходов не зави‑ сит и от номера испытания (шага), то такая последовательность называется однородной цепью Маркова с дискретным числом состояний и дискретным временем и полностью характеризу‑ ется матрицей переходов 68 5.1. Вероятностные автоматы p12 … p1k ц ч p22 … p2k ч ,    ч ч pk 2 … pkk чш где pi j — вероятность перехода из i‑го состояния в j‑е. ж p11 з p P = з 21 з  зз и pk 1 Для вероятностей p11 справедливо условие нормировки k е j =1 pi j = 1, i = 1, 2, , k . Вероятность перехода из исходного i‑го состояния в j‑е на s‑м шаге для однородной цепи вычисляется следующим образом: 1) на первом шаге Pi j (1) = pi j , 2) на втором шаге k ( ) Pi j ( 2 ) = pi 1 p1 j + pi 2 p2 j +  + pi k pkj = е pi n pnj = P 2 n =1 3) на третьем k k n =1 n =1 ( ) Pi j (3) = е Pi n ( 2 ) pnj = е P 2 in ij , ( ) pnj = P 3 ij и на s‑м шаге имеем формулу k k ( Pi j ( s ) = е Pi n ( s - 1) pnj = е P s -1 n =1 n =1 ) in ( ) pnj = P s ij . Полученный результат можно записать и в виде формулы k Pi j ( s ) = е Pi n (m ) Pnj ( s - m ) , m < n, n =1 которая называется формулой Маркова. Заметим, что в неоднородной цепи Маркова матрица пере‑ ходов зависит от номера испытания (шага). Для изучения процесса работы вероятностного автомата сле‑ Z : z1 z2  zl . Если на‑ дует задать начальное распределение p : p1 p2  pl 69 5. Дисктетно-стохастические схемы (Р‑схемы) чальное состояние zi известно, то в распределении вероятно‑ стей все нули, кроме единицы в i‑ой позиции. ПРИМЕР 23. Пусть Z = {z1 , z2 }, X = {x1 , x2 }, Y = {y1 , y2 , y3 }. Множество YP задано следующим образом: (z1, y1) (z1, y2) (z1, y3) (z2, y1) (z2, y2) (z2, y3) (x1, z1) 1 4 1 6 1 6 1 6 1 8 1 8 (x1, z2) 1 6 1 12 1 8 1 8 1 4 1 4 (x2, z1) 1 3 1 6 1 6 1 6 1 12 1 12 (x2, z2) 1 4 1 4 1 8 1 8 1 8 1 8 В начальный момент времени автомат находится в состоя‑ нии z1 . Рассмотрим процесс работы автомата при следующей последовательности входных сигналов x1 , x1 , x2 , x1 . 1‑й такт. Матрица вероятностей переходов при поступлении входно‑ го сигнала x1 имеет вид 5ц ж 7 з 12 12 ч з ч, 5ч з 3 з ч 8ш и 8 а матрица вероятностей выходов при поступлении входного сигнала x1 имеет вид ж 5 з 12 з з 7 з и 24 70 7 24 8 24 7 ц 24 чч . 9 ч ч 24 ш 5.1. Вероятностные автоматы Таким образом, система перейдет в состояния z1 и z2 с веро‑ ятностями 5ц ж 7 з 12 12 ч ж 7 5 ц ( p ( z1 ), p ( z2 )) = (1, 0 ) з 3 5 ч = з 12 , 12 ч ш з ч и з ч 8ш и 8 и сгенерирует выходные сигналы с вероятностями ж 5 з ( p ( y1 ), p ( y2 ), p ( y3 )) = (1, 0 ) з 127 з з и 24 7 24 8 24 7 ц 7 ц 24 чч = ж 5 7 . , , з 9 ч и 12 24 24 чш ч 24 ш 2‑й такт. Матрицы вероятностей переходов и вероятностей выходов такие же, как на первом такте. Тогда и 5ц ж 7 з 5 ц 12 12 ч ж 143 145 ц ч= , , ( p ( z1 ), p ( z2 )) = жз 127 , 12 чз 3 5 ч зи 288 288 чш и шз з ч 8ш и 8 ж 5 ж 7 5 ц з 12 ( p ( y1 ), p ( y2 ), p ( y3 )) = з 12 , 12 ч з 7 и шз з и 24 7 24 8 24 7 ц 24 чч = ж 105 , 89 , 96 ц . 9 ч зи 288 288 288 чш ч 24 ш 3‑й такт. Матрица вероятностей переходов теперь имеет вид ж2 з3 з з5 з и8 1ц 3 чч , 3ч ч 8ш 71 5. Дисктетно-стохастические схемы (Р‑схемы) а матрица вероятностей выходов ж1 1 1ц з2 4 4ч з ч. з3 3 1ч з ч и8 8 4ш Вероятности переходов на этом такте ж2 145 ц з 3 , ( p ( z1 ), p ( z2 )) = жз 143 чз и 288 288 ш з 5 з и8 Вероятности выходов ж1 з 143 145 ( p ( y1 ), p ( y2 ), p ( y3 )) = жз 288 , 288 цч з 32 и шз з и8 1ц 3 чч ж 4463 2449 ц = , . 3 ч зи 6912 6912 чш ч 8ш 1 4 3 8 1ц 1 ц 4 чч ж 1007 721 =з , , . 2 ч и 2304 2304 24 чш ч 8ш 4‑й такт. Матрицы вероятностей переходов и вероятностей выходов такие же, как на первом и втором тактах. Тогда и 5ц ж 7 з 4463 2449 ц 12 12 ч ж 84523 81365 ц ч= , , , ( p ( z1 ), p ( z2 )) = жз 6912 чз 5 ч зи 165888 165888 чш 6912 ш з 3 и з ч 8ш и 8 72 ж 5 з 4463 2449 ц 12 , ( p ( y1 ), p ( y2 ), p ( y3 )) = жз 6912 чз 6912 ш з 7 и з и 24 ж 20591 50833 26641 ц =з , , ч. и 55296 165888 82944 ш 7 24 8 24 7 ц 24 чч = 9 ч ч 24 ш 5.1. Вероятностные автоматы Рассмотренная вычислительная схема соответствует вычис‑ лениям в неоднородной цепи Маркова. ПРИМЕР 24. Множество YP может быть таким, что ма‑ трица переходов остается неизменной во всех тактах. Напри‑ мер, если бы в рассмотренном выше примере это множество имело вид: (z1, y1) (z1, y2) (z1, y3) (z2, y1) (z2, y2) (z2, y3) (x1, z1) 1 4 1 6 1 6 1 6 1 8 1 8 (x1, z2) 1 6 1 12 1 8 1 8 1 4 1 4 (x2, z1) 5 12 1 6 1 6 1 4 (x2, z2) 1 8 1 8 1 8 1 2 1 8 то матрица вероятностей переходов на всех тактах была бы равна 5ц ж 7 з 12 12 ч ч. Т =з 5ч з 3 з ч 8ш и 8 В этом случае вероятности Pi j ( s ) перехода из i‑го состояния в j‑е состояние на каждом такте вычисляются по формулам, исполь‑ зуемым в однородной цепи Маркова, а именно, на первом шаге Pi j (1) = Т i j , на втором шаге Pi j ( 2 ) = (Т 2 ) , и на n‑м шаге ij ( ) Pi j (n ) = T n ij . В заключение отметим, что вероятностные автоматы могут быть использованы в качестве элементов вероятностных кле‑ точных автоматов. 73 5. Дисктетно-стохастические схемы (Р‑схемы) 5.2. Задачи 1. Используя данные ПРИМЕРА 24, найти Р11(7), Р12(7), Р21(7), Р22(7). 2. Пусть Z = {z1 , z2 }, X = {x1 , x2 }, Y = {y1 , y2 }. Предложить та‑ кое множество YP , чтобы матрицы вероятностей переходов и выходов не зависели от номера такта, и найти p ( z1 ) , p ( z2 ) и p ( y1 ) , p ( y2 ) на шестом такте. 3. Рассмотреть клеточный вероятностный автомат, представ‑ ляющий собой линейную сеть из 8‑ми вероятностных автома‑ тов, каждый из которых имеет два внутренних состояния: 0 и 1. Переходы происходят по правилу: если сумма вероятностей на‑ хождения в состоянии 0 двух ближайших соседей больше или равна 0.5, то переход осуществляется с матрицей ж1 з3 Т0 = з з3 з и8 а если меньше 0.5, то с матрицей 2ц 3 чч , 5ч ч 8ш ж3 1ц з4 4ч ч. Т1 = з з4 1ч з ч и5 5ш Найти вероятность состояния 1 для элементов сети на тре‑ тьем такте. 74 6. Регрессионные модели Р егрессионные модели используются при статистиче‑ ском анализе результатов эксперимента. При этом обычно считается, что эксперимент проводится при одной или нескольких непрерывных управляемых переменных х1 , х2 , , хn (входные сигналы) и одной измеряемой случайной величины y (выходной сигнал). Регрессионные модели предла‑ гают функциональные зависимости y = f ( х1 , х2 , , хn ) , кото‑ рые выводятся из определенных принципов оптимальности и используются для прогнозирования значений y при тех зна‑ чениях х1 , х2 , , хn , которые трудно задать или при которых из‑ мерения y трудно осуществимы или вообще невозможны. Регрессионные модели можно обобщить и на случай, ког‑ да входные переменные также являются случайными перемен‑ ными. 6.1. Модель одномерной линейной регрессии Пусть в эксперименте получены пары значений ( x1 , y1 ) , ( x2 , y2 ), , ( xn , yn ), где xi — управляемая переменная, а yi — слу‑ чайная величина. Для зависимости y от x предлагается модель yi = a + bxi + ei , (6.1), 75 6. Регрессионные модели где α и β — неизвестные параметры, а yi — случайная величи‑ на — ошибка, которая обусловлена как природой самого мо‑ делируемого объекта, так и ошибкой измерения yi . При этом предполагается, что все ei статистически независимы и рас‑ пределены по одному и тому же нормальному закону N ( 0, s ) . Для определения неизвестных параметров α и β использу‑ ется метод наименьших квадратов, в котором минимизирует‑ ся величина n n f (a, b ) = е ( yi - a - bxi ) = е ( ei ) . 2 i =1 2 i =1 Из необходимых условий экстремума для функции f (a, b ) получаем систему ¶f (a, b ) ¶a ¶f (a, b ) ¶b n = е 2 ( yi - a - bxi ) ( -1) = 0 i =1 n = е 2 ( yi - a - bxi ) ( - xi ) = 0 , (6.2) i =1 или n n an + bе xi = е yi i =1 i =1 n n i =1 i =1 n a е xi + bе x = е xi yi 2 i . (6.3) i =1 Если использовать обозначения, принятые для выборочных 1 n 1 n средних x = е xi , y = е yi , то систему (6.3) можно записать n i =1 n i =1 в виде a + bx = y . ax + b x 2 = xy Решая эту систему, получаем y x 1 y b= 76 x xy x2 - x 2 , a= xy x2 x2 - x 2 = y - bx . 6.1. Модель одномерной линейной регрессии Если использовать обозначения для выборочной дисперсии D [ x ] == 1 n 1 n 2 ( xi - x ) = е xi2 - x 2 е n i =1 n i =1 и выборочной ковариации K [ x, y ] = 1 n 1 n x x y y = ( ) ( ) е i е xi yi - x Ч y , i n i =1 n i =1 то b= K [ x, y ] D [x ] . Так как α и β являются функциями случайных величин yi , то они сами тоже являются случайными величинами, и из си‑ стемы (6.3) были получены их выборочные оценки, которые да‑  и b . Полученные выше формулы для a  лее будут обозначаться a  и b можно представить в другом виде. Действительно, n n 1 n xi - x ) ( yi - y ) е ( xi - x ) yi - ( xi - x ) y е ( xi - x ) yi ( е K [ x, y ] n i =1 b = = = i =1 = i =n1 . n 1 n 2 2 2 D [x ] x x x x x x ) ( i ) ( i ) е( е е n i =1 i i =1 i =1 Если обозначить n е (x i =1 то получим n b = е (x i =1 i - x ) yi Qx i - x ) = Qx , 2 n  =y -x ,a е (x i =1 i - x ) yi Qx . Полученные выше оценки параметров регрессии обладают  и b линейно зависят от следующими свойствами: 1) a   yi , i = 1, 2, , n; 2) a и b являются несмещенными оценками па‑  ] = a и M [b ] = b; 3) дисперсии a  и b мини‑ раметров α и β т. е. M [a мальны в классе линейных несмещенных оценок. 77 6. Регрессионные модели Линейная зависимость видна из полученных формул. Несме‑ щенность оценок можно установить следующим образом. Дей‑ ствительно, n й n щ n x x y x x M y ( ) ( ) ( xi - x ) (a + bxi ) [ ] е е i ъ i i ке i i =1 i =1 i =1  й щ к ъ = M b =M = = л ы Qx Qx Qx к ъ к ъ л ы n n a е ( xi - x ) bе ( xi - x ) xi = i =1 + i =1 . Qx Qx Первое слагаемое равно нулю, а второе можно переписать в виде n bе ( xi - x ) ( xi - x ) i =1 Qx = b.  имеем Для параметра a n  щ = M [ y ] - xM йb щ = 1 (a + bx ) - x b = a. M йa е i л ы n i =1 л ы Минимальность дисперсии оценок в классе линейных несме‑ щенных оценок рассмотрена в [9].  и b Дисперсии оценок a Рассмотрим дисперсию b . Используя общую формулу для дис‑ персии линейной функции случайных величин xi , i = 1, 2, , n n n йn щ n D к е ai xi + b ъ = е ai2 D [xi ] + 2е е ai a j K йлxi , x j щы i =1 j №i л i =1 ы i =1 и учитывая, что в нашем случаеyi являются независимыми, т. е. K йл yi , y j щы = 0, i № j , получим 78 6.1. Модель одномерной линейной регрессии й n щ n 2 x x y ( ) ( xi - x ) D [ yi ] е i ъ ке i ъ = i =1 D йb щ = D к i =1 . л ы Qx Qx2 к ъ к ъ л ы Так как D [ yi ] = s 2 , то окончательно n D йb щ = л ы s 2 е ( xi - x ) i =1 Q 2 x 2 = s 2Qx s 2 = . Qx Qx2 Для дисперсии α аналогично получаем й к n  щ = D к1 y - x D йa е i л ы к n i =1 к л n щ 2 - x ) yi ъ n ж 1 x ( xi - x ) ц ж1 x2 ъ = ез чч D [ yi ] = з + з Qx Qx ъ i =1 и n и n Qx ш ъ ы е (x i =1 i ц 2 чs . ш Приводя сумму в скобках к общему знаменателю, имеем n е(x 2 i - 2 xxi + x 2 ) + nx 2 n еx 2 i  щ = ж 1 + x ц s 2 = i =1 D йa s 2 = i =1 s 2 . з ч л ы nQx nQx и n Qx ш Заметим, что β имеет нормальное распределение 2 ж s2 ц ч, N з b, з Qx чш и а параметр α — нормальное распределение n ж s 2 е xi2 з i =1 N з a, з nQx зз и ц ч ч. ч чч ш 79 6. Регрессионные модели Оценка s 2  - b x . Заметим, Для оценки s 2 используют остатки e i = yi - a i что M [ ei ] = 0, i = 1, 2, ј, n. Рассмотрим остаточную сумму ква‑ дратов n ( ) Qe = е e i i =1 2 n ( ) 2  - b x . = е yi - a i i =1 Математическое ожидание остаточной суммы квадратов вы‑ числим, предварительно преобразовав эту сумму n ( ) 2 n n Qe = е yi - y - b ( xi - x ) = е ( yi - y ) - 2b е ( yi - y ) ( xi - x ) + b i =1 2 i =1 2 n i =1 n n i =1 i =1 е (x i =1 i -x) = 2 2 2 2 = е ( yi - y ) - 2bQx b + b Qx = е ( yi - y ) - b Qx . Таким образом 2 йn щ n 2 2 M [Qe ] = M к е ( yi - y ) - b Qx ъ = е M й( yi - y ) щ - M л ы л i =1 ы i =1 Учтем, что 2 M й( yi - y ) щ = D [ yi - y ] + M 2 [ yi - y ], л ы 1 n M [ yi - y ] = a + bxi - е (a + bxi ) = b ( xi - x ) , n i =1 ж n - 1 ж 1 ц2 ц 2 D y y = s n з 2 + з1 - ч ч = s 2 (n - 1) , [ i ] е з n и n ш чш i =1 и n ( ) 2 M йb щQx = Qx D йb щ + M 2 йb щ = s 2 + Qx b2 . л ы л ы кл ъы Окончательно получаем M [Qe ] = s 2 (n - 1) + Qx b2 - s 2 - Qx b2 = s 2 (n - 2 ) . 80 йb 2 щQ . кл ъы x 2 6.1. Модель одномерной линейной регрессии Несмещенная оценка имеет вид 2 = s 2 = Qe . s n-2 Доверительные интервалы для параметров α, β и σ Статистика b - b s b - b s b - b b - b Qx = nD [ x ] Ч = Ч = 2 s s s s  s й щ D b л ы Qx имеет распределение Стьюдента с n - 2 степенями свободы T (n - 2 ) . Таким образом, с вероятностью 1 - a (здесь α — уро‑ вень значимости) имеем b - b t a (n - 2 ) < nD [ x ] < t a (n - 2 ) 1s 2 2 или b - Учитывая, что s Чt 1- a 2 (n - 2 ) nD [ x ] < b < b - t a (n - 2 ) = -t 2 1- a 2 s Ч t a (n - 2 ) 2 nD [ x ] . (n - 2), окончательно получаем доверительный интервал в следующем виде s Ч t a (n - 2 ) s Ч t a (n - 2 ) 12 b - 1- 2  . F1- a (1, n - 2 ) . s  + b x для предсказаний Использование y i = a i Пусть x произвольная точка и y = a + bx + e. В нашей модели M [ y ] = a + bx . Построим доверительный интервал для M [ y ] = a + bx . Рассмотрим дисперсию  + b x щ = D й y + b ( x - x ) щ . D й y щ = D йa л ы л ы л ы Так как y и b статистически независимы [10], то 2 D й y + b ( x - x ) щ = D [ y ] + ( x - x ) D йb щ = л ы л ы 2 2 ж ( x - x ) цч 1 2 2 s 2 1 з x x s + = s + . ( ) 2 зn Qx Qx ч i =1 n и ш n =е 83 6. Регрессионные модели Статистика a + bx - y ж 1 ( x - x )2 ц ч s з + зn Qx ч и ш имеет распределение Стьюдента T (n - 2 ), и доверительный ин‑ 2 тервал для M [ y ] = a + bx в точке х имеет вид ж 1 ( x - x )2 ц ж 1 ( x - x )2 ц  з + ч < a + bx < y + t a s з + ч. as 11зn зn Qx ч Qx ч 2 2 и ш и ш Если необходимо при данном уровне значимости построить доверительные интервалы сразу для всех значений x, то получим доверительную полосу (полоса Уоркинга — Хоттелинга) в виде y - t ж 1 ( x - x )2 ц ч, a + bx = y ± ls з + зn Qx ч и ш где l = 2Fa ( 2, n - 2 ) [9]. Проверка адекватности модели Проверка адекватности модели включает проверку всех пред‑ положений модели, а именно: 1) проверка линейности модели; 2) проверка независимости ошибок в различных точках; 3) про‑ верка нормальности законов распределения ошибок и равен‑ ства дисперсий. Проверка адекватности модели проводится на основе ана‑  - b x . лиза остатков e i = yi - y i = yi - a i 1. Качественный анализ остатков. Качественный анализ остатков проводится с использовани‑ ем графиков остатков. При этом могут получиться графики сле‑ дующих видов, представленных на рис. 6.1. 84 6.1. Модель одномерной линейной регрессии а) б) ˆε i ˆε i yi в) г) ˆε i yi ˆε i yi yi Рис. 6.1 Данные рис. 6.1, а не противоречат предположению о рас‑ пределении ошибок по закону N ( 0, s ) . Из рис. 6.1, б видно, что σ зависит от yi . Из рис. 6.1, в следует, что линейная модель y = a + bx + e неадекватна и следует использовать квадратичную модель y = a + bx + gx 2 + e. Ситуация, изображенная на рис. 6.1, г, скорее всего указывает на ошибки в вычислениях. 2. Количественный анализ остатков. Пусть в точках x1 , x2 , , xk есть повторные измерения x1 x2 … xk y11 y21 … yk 1 y12 y22 … yk 2 .     y1n1 y2n2 … yknk 85 6. Регрессионные модели ni Введем средние по группам yi = е yi j . j =1  - b x . Для этого рассмотрим Проверим близость yi к y i = a i k ( суммы Qn = е ni yi - y i i =1 ) 2 k ni и Q p = е е ( yi j - yi ) . 2 i =1 j =1 В качестве критерия адекватности выбирается статистика Qn (k - 2 ) F= , имеющая распределение Фишера F (k - 2, n - k ) . Qp (n - k ) Модель считается адекватной, если F < F1- a (k - 2, n - k ) . ПРИМЕР 25. Экспериментально установлено влияние тем‑ пературы (Х) на параметр качества продукта (Y) [10]: Х: 460, 450, 440, 430, 420, 410, 450, 440, 430. 420, 410, 400, 420, 410, 400, Y: 0.3, 0.3, 0.4, 0.4, 0.6, 0.5, 0.5, 0.6, 0.6, 0.6, 0.7, 0.6, 0.6, 0.6, 0.6. Найти параметры линейной регрессии. Определить довери‑ тельные интервалы параметров регрессии при уровне значимо‑ сти 0.1 и проверить адекватность линейной модели. Решение Для определения параметров линейной регрессии использу‑ ем процедуру linfit пакета MathCAD, которая производит под‑ гонку методом наименьших квадратов. Использование этой процедуры описано в ПРИЛОЖЕНИИ. Листинг программы. Х: = (460 450 440 430 420 410 450 440 430 420 410 400 420 410 400) Y: = (0.3 0.3 0.4 0.4 0.6 0.5 0.5 0.6 0.6 0.6 0.7 0.6 0.6 0.6 0.6) 86 X := ( 460 450 440 430 420 410 450 4406.1.430 430 410 400 420регрессии 410 400 ) Модель одномерной линейной Y := ( 0.3 0.3 0.4 0.4 0.6 0.5 0.5 0.6 0.6 0.6 0.7 0.6 0.6 0.6 0.6 ) T VX := X T VY := Y ж1ц f ( x) := з ч K := linfit ( VX , VY , f ) иxш 2.48703 ж K =з и -4.59459 ґ 10 i := 0.. 14 -3 ц ч g( t) := f ( t) ЧK ш a := 2.48703 b := -0.00459459 0.7 g( t) 0.6 0.5 VYj 0.4 0.3 0.2 400 420 440 460 480 t , VXj Рис. 6.2 Построим доверительный интервал для параметра β и про‑ верим значимость регрессии. Продолжение листинга программы. n := 15 s := stderr ( VX , VY) s = 0.09 Qx := n Чvar ( VX) = 4.933 ґ 10 Inl := b - sЧ QSt Qx 3 = -6.863 ґ 10 QSt := qt ( 0.95 , n - 2) = 1.771 -3 Inh := b + sЧ QSt Qx = -2.327 ґ 10 -3 . Таким образом, доверительный интервал имеет вид (-0.00687, - 0.00233), и гипотеза H 0 : b = 0 отвергается. Регрессия значимая. Проверка с помощью F‑критерия дает F := b 2 Qx Ч = 12.871 QFisher := qF ( 0.9 , 1 , n - 2) = 3.136 . 2 s Так как F > QFisher , то гипотеза отвергается. 87 6. Регрессионные модели Находим доверительный интервал для α. Листинг программы. n -1 (VXi)2 е aInl := a - s ЧQSt Ч i=0 n -1 е aInh := a + s ЧQSt Ч = 1.518 . n ЧQx (VXi)2 i=0 n ЧQx = 3.456 . Получаем интервал (1.518, 3.456 ) . Доверительный интервал для s 2 . sInl := s Ч ( n - 2) = 0.052 qchisq ( 0.95 , n - 2) sInh := s Ч ( n - 2) = 0.198 . qchisq ( 0.05 , n - 2) Таким образом, получаем ( 0.052, 0.198) . Анализ остатков. Resti := VYi – g(VXi) 0.2 0.1 Resti - 0.1 - 0.2 0.2 0.4 0.6 0.8 VYi Рис. 6.3 На основании графика на рис. 6.3 сделать однозначный вывод о распределении остатков по закону N ( 0, s ) затрудни‑ тельно. 88 6.1. Модель одномерной линейной регрессии Проверим адекватность с помощью F‑критерия. Распределение значений по группам имеет вид: различные значения Х: 400 410 420 430 440 450 460 значения Y: 0.6 0.6 0.6 0.6 0.4 0.3 0.3 0.6 0.7 0.6 0.4 0.6 0.5 0.5 0.6 Средние по группам: 0.6 0.6 0.6 0.5 0.5 0.4 0.3. Листинг программы. m := 7 xi := ( 400 410 420 430 440 450 460 ) T T yi := ( 0.6 0.6 0.6 0.5 0.5 0.4 0.3 ) ni := ( 2 3 3 2 2 2 1 ) ж 0.6 з з 0.6 з 0.6 з yij := 0.6 з з 0.4 з 0.3 з и 0.3 6 Qp := 0.6 T 0 ц ч 0.7 0.5 ч 0.6 0.6 ч 6 ч Qn := йnii Ч(yii - a - b Чxii)2щ = 0.019 л ы ч 0.6 0 ч i=0 0.5 0 ч ч 0 0 ш 0.4 ni i -1 е е е (yiji , j - yii)2 = 0.08 i=0 j=0 F := Qn Ч ( n - m) = 0.383 Qp Ч( m - 2) qF ( 0.9 , m - 2 , n - m) = 2.726. Так как F < Qf , то модель адекватна. Замечание. Коэффициенты простой линейной регрессии можно было вычислить, используя процедуры intercept и slope пакета MathCAD (см. ПРИЛОЖЕНИЕ). Ниже представлен ли‑ стинг программы. a := intercept ( VX , VY) = 2.487 b -3 := slope ( VX , VY) = -4.595 ґ 10 . 89 6. Регрессионные модели 6.2. Модель множественной регрессии Пусть теперь в эксперименте получены не пары значений ( x1, y1 ), ( x2 , y2 ), , ( xn , yn ), а наборы (M1, y1), (M2, y2), ..., 2 , y 2 ) , , ( M n , yn ) , где M i = ( x1i , x 2i , , xk -1i ) — точка, координаты кото‑ рой являются управляемыми переменными, а yi — случайная величина. Для зависимости y от M предлагается модель: yi = b0 + b1 x1i + b2 x2i +ј+ bk ‰1 xk -1i + ei , i = 1, 2, ј, n, (6.4) где bi i = 0, 1, ј, k - 1 — неизвестные параметры, а ei — случай‑ ные величины, относительно которых предполагается, что все ei i = 1, 2, ј, n статистически независимы и распределены по од‑ ному и тому же нормальному закону N ( 0, s ) .    В векторной форме модель имеет вид Y = X b + e, где ж y1 ц ж1 x11 x21 ј xk -11 ц ж b0 ц ж e1 ц з ч з ч з ч з ч  y 1 x12 x22 ј xk -12 ч  з b1 ч  з e2 ч Y = з 2 ч, X = з , b= , e= . з  ч з  ч з  ч з ч    зз чч зз чч зз чч зз чч и yn ш и1 x1n x2n ј xk -1n ш и bk ‰1 ш и en ш n ґ k и называется матрицей Матрица X имеет размерность  плана. Неизвестный вектор b находим методом наименьших      квадратов, т. е. минимизируем функцию f (b) = (Y - X b)T (Y - X b). Раскрывая скобки, получим     (Y - X b)T (Y - X b) =             = (Y T - bT X T )(Y - X b) = Y T Y - bT X T Y - Y T X b + bT X T X b. Каждое из слагаемых T T  Tв этой  сумме является просто числом b X Y = Y X b . и, кроме того, T       Таким образом, f (b) = Y Y - 2bT X T Y + bT X T X b. 90 6.2. Модель множественной регрессии  Рассмотрим производные слагаемых в f (b) по bi . Для второ‑ го слагаемого получаем   ¶ bT X T Y  ¶ = (b0 , b1,ј, bi ,ј, bk -1 ) X TY = ¶bi ¶bi ( ) ж y1 ц з ч  y = ( 0, 0, ј, 1, ј, 0 ) X T Y = ( xi1 , xi 2 , ј, xin ) з 2 ч . з  ч зз чч и yn ш Для третьего слагаемого имеем ж b0 ц з ч з b1 ч T T  ¶ b X Xb з ч ¶ = (b0 , b1, ј, bi , ј, bk -1 ) X T X з ч = ¶bi ¶bi з bi ч з ч зз чч и bk -1 ш ( ) ж b0 ц з ч з b1 ч з ч = ( 0, 0, ј, 1, ј, 0 ) X T X з ч+ з bi ч з ч зз чч и bk -1 ш 91 6. Регрессионные модели ж0ц з ч з0ч з ч + (b0 , b1 , ј, bi , ј, bk -1 ) X T X з ч = з1 ч з ч зз чч и0ш ж b0 ц з ч з b1 ч з ч T = ( xi1 , xi 2 , ј, xi n ) X з ч + (b0 , b1 , ј, bi , ј, bk -1 ) X b з i ч з ч зз чч и bk -1 ш ж xi1 ц з ч з xi 2 ч . з  ч з ч зx ч и in ш Так как ж b0 ц з ч з b1 ч з ч ( xi1, xi 2 , ј, xi n ) X з b ч = (b0 , b1, ј, bi , ј, bk -1 ) X T з i ч з ч зз чч и bk -1 ш ж xi1 ц з ч з xi 2 ч , з  ч з ч зx ч и in ш то окончательно в векторной форме имеем   T T  ¶ bT X T Y  b ¶ X Xb  = X TY и = 2 X T X b. ¶bi ¶bi  ¶f b   Таким образом, условие  = -2 X T Y + 2 X T X b = 0 дает систе‑ ¶b му уравнений   X T Y = X T X b, (6.6) ( ) ( () 92 ) 6.2. Модель множественной регрессии которые называются нормальными уравнениями. Если матри‑ ца невырожденная, то система (6.6) имеет решение   -1 b = ( X T X ) X T Y . Это решение дает оценку параметров модели множественной линейной регрессии методом наименьших ква‑  дратов. Далее этот вектор будем обозначать b. Случай вырожденной матрицы в данном пособии не рассма‑ тривается. Сведения об этой ситуации можно найти в [9,10].  Свойства оценок b:  1) оценки линейные по Y ,   2) оценки несмещенные, т. е. M [b] = b, 3) в классе линейных несмещенных оценок они имеют наи‑ меньшую дисперсию. Оценка для s 2 . Остаточная сумма квадратов в случае множественной ре‑ грессии имеет вид (   ОСК = Y - X b ) (Y - X b ) = Y Y - b X Y -Y  Т  Т  Т  T  T   Т X b+b XT X b =    Т  Т = Y Т Y - 2b X T Y + b X T X b.   -1 Подставляя b = X T X X T Y , получаем ( )    ОСК = Y Т Y - 2Y T X  +Y T X    = Y ТY -Y T X (( X (( X T T ) X ) X -1 ) T -1 ) T (( X T ( X ) -1 XTX XTX ) T ) -1  X TY +  X TY =     X TY = Y ТY -Y T X X T X ( ) -1  X TY =   Т    Т  = Y Т Y - b X T Y = Y Т Y - b X T X b. 93 6. Регрессионные модели Обозначим матрицу X ( X T X ) X T = A, n ґ n. Тогда математи‑ -1 ческое ожидание остаточной суммы квадратов равно n     йn щ M [ОСК ] = М йлY Т Y - Y T AY щы = М к еYi 2 - еYi Ai jY j ъ = i, j л i =1 ы n n йn щ n = М к еYi 2 (1 - Ai i ) - еYi Ai jY j ъ = е M йлYi 2 щы (1 - Ai i ) - е M [Yi ] Ai j M йлY j щы. i№ j i№ j л i =1 ы i =1 Так как  M йлYi 2 щы = D [Yi ] + M 2 [Yi ] = s 2 + M 2 [Yi ] и M [Yi ] = X b , ( ) то n n n i =1 i№ j i =1 ( ) е M йлYi 2 щы (1 - Ai i ) - е M [Yi ] Ai j M йлY j щы = е s2 + M 2 [Yi ] (1 - Ai i ) n n n n i№ j i =1 i =1 i =1 - е M [Yi ] Ai j M йлY j щы = s 2 е (1 - Ai i ) + е M 2 [Yi ] - е Ai i M 2 [Yi ] n n n n i№ j i =1 i =1 i, j - е M [Yi ] Ai j M йлY j щы = ns 2 - s 2 е Ai i + е M 2 [Yi ] - е M [Yi ] Ai j M йлY j щы. Третье и четвертое слагаемые равны. Действительно, n n    T е M [Y ] = е ( X b) = ( X b) ( X b) = b i =1 2 i n е M [Y ] A i = n  i, j е ( X b) i , j =1 94 T i   Ai j X b = X b T i i =1 и 2 ij  XT Xb M йлY j щы =   ( ) ( ) A ( X b) = ( X b) X ( X X ) T j   = bT X T X b. T T -1  XT Xb = ( ) 6.2. Модель множественной регрессии Таким образом, n n i =1 i =1 ( ( M [ОСК ] = ns 2 - s 2 е Ai i = ns 2 - s 2 е X X T X n k (X k = ns 2 - s 2 е е е i =1 j =1 l k k ( = ns 2 - е е X T X j =1 l ij (X T X -1 n ) еX jl i =1 T li ) -1 jl ) k k ) -1 ( k k ( X i j = ns 2 - s 2 е е X T X j =1 l =1 = ii X lTi = ns 2 - s 2 е е X T X j =1 l ) XT -1 n ) еX jl i =1 ) (X -1 T jl Xi j = T li X ) lj = k = ns 2 - s 2 е1 = s 2 (n - k ) . j =1 Несмещенная оценка для s 2 имеет вид 2 = ОСК = s (n - k ) (   Y -Xb )( Т   Y -Xb (n - k ) ).  Дисперсии b. Рассмотрим математическое ожидание ( )( ) ) (( й     T щ M к b - b b - b ъ. л ы   -1 Учитывая, что b = X T X X T Y и X T X ( ) -1 ) T ( = XTX ) -1 — ква‑ дратная матрица размерностью k ґ k , получаем ( )( )    T й     T щ -1 -1 M к b - b b - b ъ = M йк X T X X T Y - b X T X X T Y - b щъ = л ы л ы     -1 -1 = M й X T X X T Y - b Y T X X T X - bT щ = кл ъы (( (( ) ) (( ) )( ( ) ) ) ) 95 6. Регрессионные модели  = M [(( X T X )-1 X T YY T X ( X T X )-1     - bY T X ( X T X )-1 - ( X T X )-1 X T Y bT + bbT )] =    -1 -1 -1 = X T X X T M йлYY T щы X X T X - bM йлY T щы X X T X    -1 - X T X X T M йлY щы bT + bbT .     T Так как YY T = X b + e X b + e , то       M йлYY T щы = M йл X bbT X T + ebT X T + X beT + eeT щы = X bbT X T + s 2 E , ( ) ( ( ( ) ) )( ( ) ) где E — единичная матрица. При выводе этой формулы учте‑ но, что   1. M [ e ] = M йл eT щы = 0, 2. M йлei2 щы = D [ ei ] - M 2 [ ei ] = s 2 , ж M йe12 щ M [ e1e2 ] ј M [ e1en ] ц л ы з ч з M [e e ] M йe2 щ ј M [ e e ] ч T 2 1 2 n л 2ы ч, M йл ee щы = з з ч     з ч з й 2щ ч и M [ en e1 ] M [ en e2 ] ј M лen ы ш M йлei e j щы = M [ ei ] M йлe j щы = 0, i № j . Для остальных слагаемых имеем    T     M йлY T щы = M й X b + e щ = bT X T , M йлY щы = M й X b + e щ = X b. л ы кл ъы ( ) ( ) Тогда  -1 -1 й     T щ M к b - b b - b ъ = ( X T X ) X T X bbT X T + s 2 E X ( X T X ) л ы   -1 - b bT X T X ( X T X )     -1 - ( X T X ) X T X b bT + bbT = bbT +    -1 -1 + s 2 ( X T X ) - bbT - bbT + bbT = s 2 ( X T X ) . ( )( ) ( ( 96 ) ( ) ) 6.2. Модель множественной регрессии В результате получаем матрицу s 2 ( X T X ) , диагональные -1 элементы которой являются дисперсиями коэффициентов ре‑ грессии D йлb i щы = s 2di i , где di i — диагональный элемент матрицы (X T X ) -1 . Доверительные интервалы для коэффициентов  регрессии b Рассмотрим статистику b i - bi s di i 2 . , s= s Эта статистика имеет распределение Стьюдента с n - k сте‑ пенями свободыT (n - k ) . При уровне значимости α доверитель‑ ные интервалы для коэффициентов имеют вид b i - s di i t 1- a 2 (n - k ) < bi < b i + s di i t 1- a 2 (n - k ) . Замечание. Одновременно для всех коэффициентов bi , i = 1, 2, , k - 1 доверительные интервалы имеют место с ве‑ роятностью (1 - a ) . k Доверительный интервал для s 2 c2 1-   Т  Qe Q < s2 < 2 e , n k c ) a ( a (n - k ) 2 2   Т     где Qe = Y Т Y - b X T Y = Y Т Y - b X T X b. 97 6. Регрессионные модели Доверительный интервал для отклика Построим доверительный интервал в произвольной точке для отклика Y 0 = b0 + b1 x1 + b2 x2 +ј+ bk ‰1 xk -1 . Рассмотрим дисперсию  -1    D йY 0 щ = D й x b щ = D й x ( X T X ) X T Y щ = кл ъы кл ъы л ы T -1 -1   = x ( X T X ) X T x ( X T X ) X T s2 = (1, x1, x2 , , xk -1 ) (  = x XTX ( ) -1 ( XTX XTX ) ) -1   xT s 2 = x X T X ( Тогда ( ) -1  D йY 0 - Y 0 щ = D [Y 0 ] + D йY 0 щ = s 2 1 + x X T X л ы л ы Статистика T= Y 0 - Y 0  s2 1+ x X T X ( ( ) ( -1  xT )  xT s 2 . -1 )  xT . ) имеет распределение Стьюдента T (n - k ), и доверительный ин‑ тервал для отклика имеет вид Y 0 - t a 12 (1 + x ( X X ) x ) F1- a (1, n - k ) , то гипотеза H 0 : bi = 0 отвергается. 2. Проверка гипотезы H 0 : b1 = 0, b2 = 0, ј, bk -1 = 0. Статистика критерия выбирается в виде QR F= где n ( QR = Qy - Qe = е Yi - Y i =1 ) 2 Qe (k - 1) , (n - k )   Т   Т -Y Т Y + b X T Y = -n(Y )2 + b X T Y . Статистика имеет распределение Фишера F (k - 1, n - k ) . Если F > F1- a (k - 1, n - k ) , то гипотеза отвергается. Проверка адекватности модели Пу с т ь M 1 (1, x11 , x21 , , xk -11 ) , M 2 (1, x12 , x22 , , xk -12 ) , , M m (1, x1m , xk -12 ) , , M m (1, x1m , x2m , , xk -1m ) — различные точки. В каждой точке имеем повторные испытания 99 6. Регрессионные модели M1 Y11 Y12  Y1n1 Пусть Yi = … … …  … M2 Y 21 Y 22  Y 2n2 Mm Y m1 Ym 2 .  Y m nm m 1 ni  i = X b , n = n . Y Y , е ij е i n j =1 i =1 ( ) Статистика критерия имеет вид ( ) m ni 2 ж m ni  з е е Yi j - Y i - е е Yi j - Yi i =1 j =1 з i =1 j =1 з з F =и 2 ц ж m ni з е е Yi j - Yi ч з i =1 j =1 ч n -mч з з ч и ш ( ( (n - m ) е е ( m = (m - k ) ni i =1 j =1 Yi j - Y i m ) 2 m ni ( - е е Yi j - Yi i =1 j =1 ni е е (Y i =1 j =1 ) ij - Yi ) 2 ) 2 ц ч ч m-k ч ч ш= ) (n - m ) е n (Y -Y ) = (m - k ) е е (Y -Y ) 2 m 2 i i =1 m ni i =1 j =1 i ij i 2 . i Статистика F имеет распределение Фишера F (m - k , n - m ) . Модель считается адекватной, если F < F1- a (m - k , n - m ) . ПРИМЕР 26. Экспериментальные результаты имеют вид X i1 : 1.1, 1.5, 0.9, 1.0, 2.0, 2.1, 1.5, 1.6, 1.8, 0.7, 2.1, 2.3, X i 2 : 2.2, 0.1, 1.3, 0.5, 0.7, 3.1, 2.3, 0.4, 4.3, 0.9, 0.1, 7.6, Yi : 2.5, 4.1, 6.7, 8.1, 4.0, 1.7, 1.1, 1.5, 3.9, 6.1, 7.0, 7.2 . 100 6.2. Модель множественной регрессии Найти коэффициенты модели Yi = b0 + b1 X i1 + b2 X i 2 + ei , по‑ строить доверительные интервалы для параметров регрессии, проверить гипотезы о коэффициентах регрессии и значимость модели. Листинг программы. Y := ( 2.5 4.1 6.7 8.1 4.0 1.7 1.1 1.5 3.9 6.1 7.0 7.2 ) T T ж 1 1 1 1 1 1 1 1 1 1 1 1 ц X := з 1.1 1.5 0.9 1.0 2.0 2.1 1.5 1.6 1.8 0.7 2.1 2.3 ч з ч и 2.2 0.1 1.3 0.5 0.7 3.1 2.3 0.4 4.3 0.9 0.1 7.6 ш n := 12 k := 3 b ( T := X X ) -1 ж 5.819 ц X Y = з -1.05 ч з ч и 0.153 ш T Таким образом, b0 = 5.810, b1 = -1.05, b0 = 0.153. Найдем оценку для s 2 и дисперсии коэффициентов. T T ( T Qe := Y Y - Y X X X ( T d := X X ) -1 )- 1 XT Y = 64.303 s := Qe = 2.673 n-k ж 0.884 -0.546 0.023 ц з ч = -0.546 0.409 -0.045 з ч и 0.023 -0.045 0.024 ш 2 2 Db0 := s Чd0 , 0 Db0 = 6.316 Db1 := s Чd1 , 1 Db1 = 2.924 2 Db2 := s Чd2 , 2 Db2 = 0.169 Видно, что дисперсии имеют большие значения. Находим границы доверительных интервалов. 101 6. Регрессионные модели b0l := b0 - qt ( 0.95, n - k ) Чs Ч d0 , 0 = 1.212 b0h := b0 + qt ( 0.95, n - k ) Чs Ч d0 , 0 = 10.426 b1h := b1 + qt ( 0.95, n - k ) Чs Ч d1 , 1 = 2.085 b1l b2h b2l := b1 := := - qt ( 0.95, n - k ) Чs Ч d1 , 1 = -4.184 + qt ( 0.95, n - k ) Чs Ч d2 , 2 = 0.907 b2 b2 - qt ( 0.95, n - k ) Чs Ч d2 , 2 = -0.601 Из интервала для b0 (1.212, 10.426 ) следует, что гипотеза H 0 : b0 = 0 отвергается. Из интервалов для b1 ( -4.184, 2.085) и b2 ( -0.601, 0.907 ) следует, что гипотезы H 0 : b1 = 0 и H 0 : b2 = 0 не отвергаются. Проверим значимость модели, т. е. гипотезу H 0 : b1 = 0, b2 = 0. n -1 Y i Ym := е i=0 n 2 T T = 4.492 Qr := -n ЧYm + b ЧX ЧY = 2.766 Fq := qF ( 0.9 , k - 1 , n - k ) = 3.006 F := Qr 2 s ( k - 1) = 0.194 . Так как значение критерия 0.194 меньше 3.285, то гипотеза не отвергается и регрессия незначимая. 6.3. Одномерные полиномиальные модели 1. Квадратичная регрессия Рассмотрим зависимость y = a + bx + gx 2 + e. Для эксперимен‑ тальных точек ( xi , yi ) , i = 1, 2, , n математическая модель ква‑ 102 6.3. Одномерные полиномиальные модели дратичной регрессии имеет вид y i = a + bxi + gxi2 + ei , где ei — слу‑ чайные величины, распределенные по закону N (0, s). Коэффициенты a, b, g находим по методу наименьших квадра‑ n тов, минимизируя функцию f (a, b, g ) = е ( yi - a - bxi - gxi2 ) . 2 i =1 В результате получаем систему n n n ¶f = е yi - na - bе xi - g е xi2 = 0 ¶a i =1 i =1 i =1 n n n n ¶f = е xi yi - a е xi - bе xi2 - g е xi3 = 0 . ¶b i =1 i =1 i =1 i =1 n n n n ¶f = е xi2 yi - a е xi2 - bе xi3 - g е xi4 = 0 ¶g i =1 i =1 i =1 i =1 Решение запишем в виде ж з n з з n где A = з е xi з i =1 з n 2 з е xi и i =1 n еx i i =1 n еx i =1 n 2 i еx i =1 3 i ж n ц з е yi ч з i =1 ч жaц з n ч з ч -1 з b ч = A з е xi yi ч , з i =1 ч зgч и ш з n 2 ч з е xi yi ч и i =1 ш n ц xi2 ч е i =1 ч n ч 3 xi ч . е i =1 ч n ч xi4 ч е i =1 ш Этот же результат, естественно, получается из формул мно‑ жественной регрессии при выборе матрицы плана в виде 103 6. Регрессионные модели ж1 з 1 X =з з зз и1 Тогда ж з n з з n T X X = з е xi з i =1 з n 2 з е xi и i =1 x1 x12 ц ч x2 x22 ч .   ч ч xn xn2 чш ц ч i =1 i =1 ч n n ч 2 3 xi е xi ч , е i =1 i =1 ч n n ч xi3 е xi4 ч е i =1 i =1 ш ж n ц еy ж y1 ц зз i =1 i чч ж1 1 ј 1 ц з ч  з з n ч ч y T X Y = з x1 x2 ј xn ч з 2 ч = з е xi yi ч з ч ч з x 2 x 2 ј x 2 ч з  ч з i =1 n шз 2 и 1 ч з n y и n ш з x 2 y чч е i i и i =1 ш -1 n n n ц ж 2ц ж n x x е е i i ч з е yi ч з i =1 i =1 з ч з i =1 ч ж b0 ц n n n n  з ч з ч -1 з ч и з b1 ч = X T X X T Y = з е xi е xi2 е xi3 ч з е xi yi ч . i =1 i =1 з i =1 ч з i =1 ч зb ч и ш з n 2 n 3 n 4ч з n 2 ч з е xi е xi е xi ч з е xi yi ч i =1 i =1 и i =1 ш и i =1 ш ( n n еx еx i 2 i ) ПРИМЕР 27. Используем данные ПРИМЕРА 25 для ква‑ дратичной модели. Применим формулы множественной ре‑ грессии. Листинг программы является продолжением листин‑ га ПРИМЕРА 25. Транспонированная матрица плана квадратичной модели обозначена Xsq. 104 6.3. Одномерные полиномиальные модели Листинг программы. i := 0 .. 14 Xsq0 , i := 1 Xsq1 , i := VXi Xsq2 , i := VXi ЧVXi ж 2.281 ґ 104 ц 0.124 -106.672 з ч -1 T -4ч D := ( Xsq ЧXsq ) D = з -106.672 -5.824 ґ 10 ч 0.499 з з 0.124 -4 -7 ч 6.801 ґ 10 -5.824 ґ 10 и ш b ж -17.266 ц з ч 0.088 := D ЧXsq ЧVY = з ч Gs ( t) := з -4ч и -1.079 ґ 10 ш 14 Qe := е T ж1ц з ч з t ч Чb з 2ч иt ш (VYi - Gs (VXi))2 = 0.088 i=0 ssq := Qe = 0.086 . n-3 На рис. 6.4 представлены графики экспериментальных значе‑ ний, линия линейной регрессии и линия квадратичной регрессии. 0.8 Gs ( t) VYj 0.6 0.4 g( t) 0.2 400 420 440 460 480 t , VXj , t Рис. 6.4 105 6. Регрессионные модели На рис. 6.5 представлен график остатков j := 0 .. 14 Restsqj := VYj - Gs ( VXj ) 0.2 0.1 Restsqj - 0.1 - 0.2 0.2 0.3 0.4 0.5 0.6 0.7 VYj Рис. 6.5 Доверительные интервалы для полученных коэффициентов регрессии и проверку гипотез рассмотрим далее после исследо‑ вания общей полиномиальной модели. 2. Общая одномерная полиномиальная регрессия Полиномиальная модель в общем виде y = b0 + b1 x + b2 x 2 +ј+ bk -1 x k -1 + b1 x + b2 x 2 +ј+ bk -1 x k -1 + e содержит k неизвестных коэффициентов. В точ‑ ках x1 , x2 , , xn имеем yi = b0 + b1 xi + b2 xi2 +ј+ bk -1 xik -1 + ei , i = 1, 2, ј, k -1 + ei , i = 1, 2, ј, n, n > k . Используя матрицу плана k -1 xi ж1 з 1 X =з з зз и1 x1 x12 x2 x22   xn xn2 ј ј ј ј x1k -1 ц ч x2k -1 ч , ч  ч xnk -1 чш решение получаем в виде  -1  b = ( X T X ) X T y. 106 6.3. Одномерные полиномиальные модели Замечание. С увеличением порядка регрессионной модели (с ростом k) матрица плана становится плохо обусловленной и возрастает влияние ошибок округления, в силу чего можно получить совершенно неверный результат. Доверительные интервалы Доверительные интервалы вычисляются по тем же форму‑ лам, что и для множественной регрессии b i - s di i t 1- a 2 (n - k ) < bi < b i + s di i t 1- a 2 (n - k ) , Qe Q < s2 < 2 e . c a (n - k ) c a (n - k ) 2 1- 2 2 Проверка значимости коэффициентов при старших членах Пусть есть две модели. Модель 1: yi = b0 + b1 xi + b2 xi2 +ј+ bk -1 xik -1 + ei . Модель 2: yi = b0 + b1 xi + b2 xi2 +ј+ bq -1 xiq -1 + ei , q < k . Проверим гипотезу H 0 : bq = 0, bq +1 = 0, ј, bk -1 = 0. Пусть S1 и S 2 — остаточные суммы первой и второй модели, соответ‑ ственно. ( n - k ) ( S 2 - S1 ) , которая имеет Рассмотрим статистику F = ( k - q ) S1 распределение Фишера F (k - q, n - k ) . Если F > F1- a (k - q, n - k ) , то гипотеза отклоняется. 107 6. Регрессионные модели Пример 27 (продолжение) Рассмотрим значимость коэффициента b2 . В этом случае до‑ верительный интервал имеет вид b 2 - t0,95 (n - 3) s d22 < b2 < b 2 + t0,95 (n - 3) s d22 . Приведем результаты расчета границ доверительного ин‑ тервала. Продолжение листинга. b2 + ( qt ( 0.95 , n - 3) ) Чssq Ч D2 , 2 = 1.798 ґ 10 b2 - qt ( 0.95 , n - 3) Чssq Ч D2 , 2 = -2.338 ґ 10 -5 -4 Таким образом, гипотеза H 0 : b2 = 0 не отвергается. С помо‑ щью F‑критерия получаем F := ( n - 3) Ч ( s - ssq) = 0.572 , ( 3 - 2) Ч s qF ( 0.9 , 1 , n - 3) = 3.177 что подтверждает полученный выше результат. Проверим адекватность квадратичной модели. Вычисляем остаточные суммы. Изменяется только остаточная сумма Qn. Ниже представлены результаты вычисления. 6 Qnsq := е й(yii - Gs (xii))2 Чniiщ = 5.423 ґ 10- 3 л ы i=0 6 Qp := ni i -1 е е (yiji , j - yii)2 = 0.08 i=0 j=0 F := Qnsq Ч ( 15 - 7) = 0.136 qF ( 0.9 , 4 , 8) = 2.806 . Qp Ч( 7 - 3) Таким образом, квадратичная модель также оказывается адекватной. 108 6.4. Планирование регрессионного эксперимента Рассмотрим общую схему выбора порядка полиномиаль‑ ной модели. 1. Сравниваем остаточные суммы квадратов. Идеально, если с ростом k остаточная сумма квадратов уменьшается и выходит на некоторое постоянное значение. 2. Проверяем значимость коэффициентов при старших чле‑ нах. Если обнаружен первый незначимый старший коэффици‑ ент, то увеличение степени прекращаем. При этом следует учи‑ тывать четность или нечетность функции регрессии. 3. Используем обратный отбор, т. е. задаем полином с мак‑ симальным значением k и, последовательно уменьшая это зна‑ чение, проверяем адекватность. Правда, в этом случае неясно, с какого именно максимального значения следует начать эту процедуру, т. е. фактически обратный отбор проходит в предпо‑ ложении верности предположения о максимальном значении k. 6.4. Планирование регрессионного эксперимента Общая стратегия планирования регрессионного экспери‑ мента включает: 1) выбор области изменения переменных, 2) выбор расположения точек в этой области, 3) выбор регрессионной модели, 4) определение числа повторных наблюдений и точек, в ко‑ торых следует проводить повторные наблюдения. Выбор области изменения переменных определяется, как правило, условиями проведения эксперимента. Выбор распо‑ ложения точек зависит от цели, которую преследует экспери‑ ментатор. Например, если нужно минимизировать дисперсии  коэффициентов регрессии D[b ], то используют D-оптимальные планы. 109 6. Регрессионные модели Рассмотрим D‑оптимальный план для модели y = a + bx + e. Пусть область изменения переменной xi О [ x min , x max ]. Изменим x + x max xi - min 2 масштаб так, чтобы все точки xi = принадлежали x min - x max 2 отрезку [ -1, 1]. Тогда минимальное значение дисперсии углово‑ s2 го коэффициента D йлb щы = n будет достигаться при мак‑ 2 x x ( ) е i i =1 n симальном значении Qx = е ( xi - x ) . Максимальное значение 2 i =1 этой величины будет в случае, когда точки поровну поделены между концами отрезка. Замечание. Несмотря на наличие повторных измерений, проверить адекватность модели не удастся, так как число ко‑ эффициентов регрессии равно числу различных точек. Таким образом, этот план следует использовать только, если есть пол‑ ная уверенность в выбранной модели. В случае общей одномерной полиномиальной модели y = b9 + b1 x + b2 x 2 +ј+ bk -1 x k -1 + e D‑оптимальное распределение содержит по npi точек при каждом xi , где n — общее число то‑ 1 чек, а pi = [9]. Точки xi находятся из решения уравнения k n -1) d n ( 2 (1 - x ) Pkў-1 ( x ) = 0, где Pn ( x ) = 2n n! dx n (1 - x 2 ) — полином Лежан‑ дра n‑го порядка. Замечание. Если число npi окажется дробным, то проводит‑ ся округление. 110 6.4. Планирование регрессионного эксперимента ПРИМЕР 28. Рассмотрим модель y = b9 + b1 x + b2 x 2 + e. В этом ( -1) d 2 1 1 1 - x 2 = -1 + 2 x 2 . Урав‑ случае k = 3, pi = и P2 ( x ) = 2 2 2 3 2 2! dx 2 1 x 6 x = , x нение имеет вид а его корни 1,2 = ±1, x3 = 0. Пусть 2 ( ( ) ( ) ) общее число точек равно 6. Тогда получаем по две точки при каждом значении xi . Заметим, что оптимальный план для од‑ ной модели становится неоптимальным для другой. ПРИМЕР 29. Рассмотрим модель y = b9 + b1 x + b2 x 2 + b2 x 3 + e. 1 1 Имеем k = 4, pi = и P3 ( x ) = ( 9 x - 5x 3 ) . Уравнение имеет вид 4 2 (1 - x 2 ) (9 - 15x 2 ) = 0, а его корни x1,2 = ±1, x3,4 = ± 35 . Если целью является минимизировать дисперсию только старшего коэффициента, т. е. D йb k -1 щ ® min, то оптимальный л ы план предусматривает выбор точек по формуле ж ip ц xi = - cos з ч , i = 0, 1, ј, k - 1. и k -1 ш Вероятности, с которыми наблюдения распределяются по этим точкам, находятся по формуле 1 м п 2 (k - 1) , i = 0, i = k - 1 п pi = н п 1 , i = 1, 2, , k - 2 по k - 1 Если целью является минимизировать дисперсию прогноза й D b 0 + b 1 x + b 2 x 2 + … + b k -1 x k -1 щ , x > 1, то план имеет вид л ы Li ( x ) ж ip ц xi = - cos з , i = 0, 1, , k - 1, pi = k -1 , ч и k -1 ш L x е i( ) i =0 111 6. Регрессионные модели где Li ( x ) = ( x - x0 ) ( x - x1 ) ( x - xk -1 ) . ( xi - x0 ) ( xi - x1 ) ( xi - xi -1 ) ( xi - xi +1 ) ( xi - xk -1 ) ПРИМЕР 30. Рассмотрим оптимальный план для прогноза простой линейной регрессии. В этом случае i = 0 и i = 1, ж 0Чp ц ж 1Ч p ц = -1, x1 = - cos з x0 = - cos з ч ч = +1. Для вероятностей по‑ и 1 ш и 1 ш лучаем значения L0 ( x ) = p0 = p1 = L0 L0 + L1 L1 ( x - x1 ) , ( x0 - x1 ) = = L1 ( x ) = x - x1 x0 - x1 x - x0 x - x1 + x0 - x1 x0 - x1 x - x0 x0 - x1 = ( x - x0 ) , ( x1 - x0 ) x - x1 x - x1 + x - x0 = x - x0 , . x - x1 + x - x0 x - x0 x - x1 + x0 - x1 x0 - x1 0.2 1 2.2 11 = , p1 = = Пусть х = 1, 2 и n = 12. Тогда p0 = 0.2 + 2.2 12 0.2 + 2.2 12 и оптимальный план содержит одну точку при х = -1 и один‑ надцать точек при х = +1. Для квадратичной регрессии имеем L0 + L1 ж 0Чp ц i = 0 и x0 = - cos з ч = -1, и 2 ш ж 1Ч p ц i = 1 и x1 = - cos з ч = 0, и 2 ш 112 6.4. Планирование регрессионного эксперимента ( x - x1 ) ( x - x2 ) ж 2Чp ц = +1. L0 ( x ) = i = 2 и x2 = - cos з , ч ( x0 - x1 ) ( x0 - x2 ) и 2 ш L1 ( x ) = ( x - x0 ) ( x - x2 ) , ( x1 - x0 ) ( x1 - x2 ) L2 ( x ) = ( x - x0 ) ( x - x1 ) . ( x2 - x0 ) ( x2 - x1 ) Для вероятностей получаем p0 = p1 = p2 = L0 L0 + L1 + L2 L1 L0 + L1 + L2 L2 L0 + L1 + L2 = ( x - x1 ) ( x - x2 ) ( x0 - x1 ) ( x0 - x2 ) , ( x - x1 ) ( x - x2 ) + ( x - x0 ) ( x - x2 ) + ( x - x0 ) ( x - x1 ) ( x0 - x1 ) ( x0 - x2 ) ( x1 - x0 ) ( x1 - x2 ) ( x2 - x0 ) ( x2 - x1 ) = ( x - x0 ) ( x - x2 ) ( x1 - x0 ) ( x1 - x2 ) , ( x - x0 ) ( x - x2 ) + ( x - x0 ) ( x - x1 ) ( x1 - x0 ) ( x1 - x2 ) ( x2 - x0 ) ( x2 - x1 ) = ( x - x1 ) ( x - x2 ) + ( x0 - x1 ) ( x0 - x2 ) ( x - x0 ) ( x - x1 ) ( x2 - x1 ) ( x2 - x0 ) . ( x - x1 ) ( x - x2 ) + ( x - x0 ) ( x - x2 ) + ( x - x0 ) ( x - x1 ) ( x0 - x1 ) ( x0 - x2 ) ( x1 - x0 ) ( x1 - x2 ) ( x2 - x0 ) ( x2 - x1 ) Для х = 1, 2 получаем L0 = 0.12, L1 = 0.44, L2 = 1.32 и p0 = 0.12 3 0.44 11 1.32 33 = , p2 = = . = , p1 = 1.88 47 1.88 47 1.88 47 113 6. Регрессионные модели Множественная линейная регрессия Для модели множественной линейной регрессии Yi = b0 + b1 x1i + b2 x2i +ј+ bk -1 xk -1i + ei c матрицей ж1 x11 … xk -11 ц ч з 1 x12 … xk -12 ч X =з з    ч чч зз и1 x1n … xk -1n ш справедлива теорема: n n  пусть е xi j = 0 и е xi2j = C , тогда минимум D[b] достигается, j =1 j =1 если столбцы матрицы X взаимно ортогональны [9]. 6.5. Задачи 1. Защитное покрытие деталей содержит некоторое количе‑ ство железа в процентах (Х) и испытывается на коррозию в те‑ чение определенного времени в агрессивной среде. Коррозия определяется по количеству разрушенного покрытия на ква‑ дратный сантиметр в процентах (Y). Х: 0.1, 0.2, 0.15, 0.12, 0.3, 0.26, 0.33, 0.18, 0.1, 0.3, 0.2, 0.15 Y: 0.3, 0.5, 0.5, 0.28, 0.7, 0.7, 0.8, 0.52, 0.2, 0.6, 0.45, 0.4 Найти параметры простой линейной модели. Определить их доверительные интервалы при уровне значимости α = 0.05, проверить значимость регрессии. Проверить адекватность ква‑ дратичной модели при том же уровне значимости. 2. Стоимость ремонта автомобиля (Y тыс. руб.) в год в за‑ висимости от возраста автомобиля (X лет) представлена ниже. X: 1.0, 1.5, 2.0, 3.0, 4.5, 5.0, 6.0, 7.0, 1.5, 3.0, 5.0, 6.0, 4.5 Y: 5, 3, 14, 21, 40, 55, 70, 60, 5, 25, 50, 60, 47 114 6.5. Задачи Найти коэффициенты линейной и квадратичной регрессии, их доверительные интервалы при α = 0.05 и проверить адекват‑ ность моделей. 3. Экспериментатор наблюдает изменение линейного разме‑ ра устройства в зависимости от температуры. Получены следу‑ ющие результаты: T: 50, 75, 100, 125, 150, 175, 200, 160, 80, 140, 60, 55 D: 4.5, 5.0, 5.0, 5.2, 5.8, 6.5, 8.0, 7.0, 4.8, 5.5, 4.0, 4.0 Далее он получил оценки параметров линейной модели b0 , b1 , s 2 , прогноза в точке T = 220, вычислил их дисперсии, про‑ верил значимость и адекватность модели. Потом по какой-то причине он решил, что модель должна быть только линейной, и получил оценки всех параметров из оптимального плана: Т = 50 Е = 200 4.5 8.1 4.3 7.9 3.9 8.2 4.4 8.0 3.8 8.4 4.6 7.9 Сравнить результаты обоих планов. 4. Экспериментатор наблюдает изменение линейного разме‑ ра устройства в зависимости от температуры. Получены следу‑ ющие результаты: T: 50, 75, 100, 125, 150, 175, 200, 160, 80, 140, 60, 55 D: 4.5, 5.0, 5.0, 5.2, 5.8, 6.5, 8.0, 7.0, 4.8, 5.5, 4.0, 4.0 Далее он получил оценки параметров квадратичной модели b0 , b1 , s 2 , прогноза в точке T = 220, вычислил их дисперсии, прове‑ рил значимость и адекватность модели при α = 0.05. Потом по какойто причине он решил, что модель действительно квадратичная, и по‑ лучил оценки всех параметров из оптимального плана: Т = 50 Е = 125 Т = 200 4.5 5.1 8.1 4.3 5.5 7.9 3.9 5.6 8.2 4.4 4.9 8.0 Сравнить результаты обоих планов. 115 6. Регрессионные модели 5. Исследуется влияние температуры (Т) и давления (Р) на время работы прибора до отказа (t). Получены результаты: T P t 80 0.8 100 120 1.0 100 150 0.9 80 200 1.0 70 120 3.0 50 150 2.0 60 150 0.9 90 200 1/0 65 80 0.8 130 Найти коэффициенты множественной линейной регрессии. Определить их доверительные интервалы при уровне значимо‑ сти a = 0.05, проверить гипотезы о параметрах. 116 7. Непрерывные стохастические модели (Q‑схемы) 7.1. Системы массового обслуживания Q ‑схемы являются математическими схемами теории массового обслуживания. Процессы в системах мас‑ сового обслуживания (СМО) могут иметь различную физическую природу. Это могут быть экономические системы, производственные системы, технические системы, информа‑ ционные и компьютерные системы. Характерными особенно‑ стями СМО являются: 1) случайный характер поступления заявок (требований) на обслуживание; 2) случайный характер длительности обслуживания. Работа СМО включает два необходимых элемента: а) ожида‑ ние обслуживания и б) собственно обслуживание, которое пред‑ полагает наличие прибора обслуживания. На рис. 7.1 приведе‑ на элементарная схема работы одноканальной СМО. Y U К Н X Рис. 7.1 117 7. Непрерывные стохастические модели (Q‑схемы) Здесь Х — входной поток, Н — накопитель с заданными па‑ раметрами (например, с заданной емкостьюCn), U — поток об‑ служивания, К — канал обслуживания, Y — выходной поток. Таким образом, модель СМО есть совокупность Q = X , H ,U , Z , A , где Z = {z H , z K }, z H — число заявок в накопи‑ теле, z K — число заявок в канале обслуживания. Обычно z K = 0, если в канале нет заявок, и z K = 1, если заявки есть. А — алго‑ ритмы обслуживания. Различные типы СМО отличаются друг от друга по следую‑ щим параметрам: 1) по типу ожидания: 1.1) время ожидания не ограничено, 1.2) время ожидания ограничено, 1.3) СМО без ожидания (если канал закрыт, то заявка не принимается); 2) по типу постановки в очередь: 2.1) по очереди в порядке поступления, 2.2) с учетом приоритета; 3) по типу каналов обслуживания: 3.1) одноканальные СМО, 3.2) многоканальные СМО; 4) по длительности обслуживания: 4.1) с неограниченным временем обслуживания, 4.2) с ограниченным временем обслуживания. Рассмотрим моделирование входного потока. Пусть началь‑ ный момент времени t0 = 0, и в случайные моменты {t1 ,t2 ,,tn ,} поступают заявки. Заметим, что промежуток между поступлени‑ ями заявок Dti = ti +1 - ti , естественно, тоже случайная величина. Поток называется однородным, если поступают заявки од‑ ного типа. Такой поток полностью характеризуется набором {t0 ,t1,t2 ,,tn ,}. Если поступают заявки с различным набором признаков, то такой поток называется неоднородным. 118 7.2. Марковские процессы Поток называется ординарным, если в каждый момент вре‑ мени ti поступает только одна заявка, или, если p — вероятность поступления двух или более заявок за промежуток Dt , то p ® 0 при Dt ® 0. Поток называется стационарным, если число заявок, посту‑ пивших за промежуток Dt , зависит только от длины Dt и не за‑ висит от расположения этого отрезка на оси. 7.2. Марковские процессы В рассмотренных в разделе 5 марковских цепях время име‑ ло дискретный характер. Однако в марковской цепи время мо‑ жет быть непрерывным. В этом случае имеем марковский по‑ ток x (t ) , где x (t ) описывает состояние системы в момент t. Рассмотрим характеристики этого процесса. Пусть в момент s система находится в состоянии i и pi j ( s,t ) — вероятность того, что в момент t > s она будет в j‑м состоянии, независимо от того, в каком состоянии она была до момента s. Если pi j ( s,t ) = pi j (t - s ) , то марковский процесс называется однородным. Переходные вероятности pi j ( s,t ) удовлетворяют уравнению Колмогорова k pi j ( s,t ) = е pi n ( s,u )pn j (u,t ) , s Ј u Ј t . n =1 В случае однородного процесса k pi j (t - s ) = е pi n (u - s )pn j (t - u ) , n =1 и с точностью до обозначений имеем уравнение Колмогоро‑ ва — Чепмена k pi j ( s + t ) = е pi n ( s )pn j (t ) , n =1 i, j = 1, 2, ..., k. 119 7. Непрерывные стохастические модели (Q‑схемы) Если переходные состояния обладают свойством непре‑ м1, i = j pi j (t ) = н рывности pi j ( 0 ) = lim и непрерывно дифферен‑ t ®0 о0, i № j цируемы при t > 0, то существуют конечные пределы lim pi j (t ) - pi j ( 0 ) = li j , где li j есть плотности перехода из состоя‑ t t ®0 ния i в состояние j. Если система не изменяет свое состояние бо‑ лее одного раза за произвольно малый промежуток времени, то k еl j =1 k ij = 0 и е li j = -li i = li [11]. В этом случае марковский про‑ j №i цесс полностью характеризуется матрицей плотности li j , i, j = 1, 2, ј, k и, если в момент t = 0 система находилась li j , i № j есть вероятность того, что при вы‑ в состоянии i, то p = li ходе из i‑го состояния система перейдет именно в j‑е состояние. При сделанных выше предположениях переходные вероят‑ ности удовлетворяют следующим системам дифференциаль‑ ных уравнений: 1. Пусть P (t ) = pi j (t ) — матрица переходных состояний. Тогда P ў (t ) = lim P (t + Dt ) - P (t ) Dt Dt ® 0 = P (t ) lim = lim Dt ® 0 P ( Dt ) - E Dt ® 0 Dt P (t ) P ( Dt ) - P (t ) Dt = = P (t ) L, где L = li j , i, j = 1, 2, ј, k . Система P ў (t ) = P (t ) L называется прямой системой уравне‑ ний Колмогорова и имеет вид k piў j = е pi n l n j , i, j = 1, 2, ј, k . n =1 120 7.2. Марковские процессы Из системы P ў (t ) = P (t ) L следует P( и Ґ P (t ) - P ( 0 ) = е P ( m) m) (t ) = P (t ) Lm (0 ) m =1 Ґ L mt m Ґ L mt m tm = еE =е . m ! m =1 m! m =1 m ! 2. Если при вычислении производной представить предел в виде P ( Dt ) P (t ) - P (t ) P (t + Dt ) - P (t ) P ў (t ) = lim = lim = Dt ® 0 Dt ® 0 Dt Dt ( P ( Dt ) - E ) P (t ) = LP t , = lim () Dt ® 0 Dt то получим обратную систему уравнений Колмогорова k P ў (t ) = LP (t ) или piў j = е li n pnj , i, j = 1, 2, , k . n =1 Обе системы имеют одинаковые решения. Пуассоновский процесс Пуассоновский процесс является примером марковского процесса, используемым для моделирования входного потока в системах массового обслуживания. Пусть x (t ) — число собы‑ тий, наступивших в момент t. Пусть x (t ) — поток независимых событий, обладающий следующими свойствами: 1) вероятность отдельного события за время Dt равна lDt + o ( Dt ), где λ — положительная постоянная; 2) вероятность наступления более, чем одного события за время Dt есть o ( Dt ); 121 7. Непрерывные стохастические модели (Q‑схемы) 3) если Dt1 , Dt2 , , Dtn — непересекающиеся временные ин‑ тервалы, то количества событий x ( Dt1 ) , x ( Dt2 ) , ј, x ( Dtn ) есть взаимно независимые случайные величины. При этих условиях пуассоновский процесс является одно‑ мl, j = i + 1 родным марковским процессом и li j = н . о0, j № i + 1 Заметим, что число состояний системы не ограничено, ма‑ трица плотностей перехода является треугольной ж -l з з0 L = з0 з зј з0 и l -l ј l -l ј ј ј ј  ј 0 ц ч 0 ч 0 ч ч l ч -l чш и система уравнений Колмогорова имеет вид: Ґ piў j (t ) = е li n pnj (t ) = -li i pi j (t ) + li i +1 pi +1 j (t ) . n=0 Учитывая, что pi j(t) = p0 j–i(t), получаем Ґ p0ў j -i (t ) = е l 0 n pn j -i (t ) = l 0 0 p0 j -i (t ) + l 01 p1 j -i (t ) = n=0 = -lp0 j -i (t ) + lp1 j -i (t ) . Обозначим p0 j -i (t ) = Pj (t ) , тогда система преобразуется к виду мпP0ў (t ) = -lP0 (t ) , j = 0 н опPjў (t ) = -lPj (t ) + lPj -1 (t ) , j = 1, 2, ј , с начальными условиями P0 ( 0 ) = 1, Pj ( 0 ) = 0, j = 1, 2,  . Система имеет решение Pj (t ) = (lt ) j! j e - lt , что совпадает с вероятностью на‑ ступления j событий в распределении Пуассона с параметром lt. 122 7.3. Модель одноканальной СМО с отказами 7.3. Модель одноканальной СМО с отказами ПРИМЕР 31. Одноканальная система массового обслужи‑ вания с отказами. Математическая модель Пусть одноканальная система массового обслуживания мо‑ жет находиться в двух состояниях: 0 — система свободна; 1 — система занята. На обслуживание поступает пуассоновский по‑ ток заявок с параметром λ. Время обслуживания есть случайная величина τ, распределенная по экспоненциальному закону t P ( t < t ) = т me - mx dx = 1 -e - mt , с параметром μ. Пусть в момент t1 си‑ стема была занята и τ — время окончания обслуживания. Тог‑ - m (t -t1 ) , и за промежуток t - t1 система да при t > t1 имеем P ( t > t ) = e - m (t -t1 ) переходит в состояние 0. с вероятностью 1 - e Переходные вероятности связаны соотношениями p00 = 1 - p01 , p11 = 1 - p10 , и обратная система уравнений Колмо‑ горова приобретает вид ў (t ) = p00 (t ) l 00 + p01 (t ) l10 = p00 = - p00 (t ) l + (1 - p00 (t )) m = - ( l + m ) p00 (t ) + m, ў (t ) = p10 (t ) l 01 + p11 (t ) l11 = p11 = (1 - p11 (t )) l - p11 (t ) m = - ( l + m ) p11 (t ) + l. Видно, что система распадается на отдельные линейные неоднородные уравнения, решения которых p00 (t ) = С1e - ( l + m )t + m l - l +m t , p11 (t ) = С 2e ( ) + . l +m l +m 123 7. Непрерывные стохастические модели (Q‑схемы) Из начальных условий p00 ( 0 ) = 1, p11 ( 0 ) = 1 следует С1 = 1 - m l , С2 = 1 . l +m l +m Таким образом, p00 (t ) = l - ( l + m )t m m - ( l + m )t l e + , p11 (t ) = e + . l +m l +m l +m l +m При t ® Ґ вероятности экспоненциально стремятся к стацио­ m l , p11 = . Заметим, что стацио‑ нарным значениям p00 = l +m l +m нарные значения вероятностей (если они существуют) удовлет‑ воряют системе P L = 0. Для получения решения можно было использовать и пря‑ мую систему уравнений Колмогорова, которая в данном слу‑ чае имеет вид: ў (t ) = l 00 p00 (t ) + l 01 p10 (t ) = p00 = -lp00 (t ) + l (1 - p11 (t )) - = l - l ( p00 (t ) + p11 (t )) , ў (t ) = l10 p01 (t ) + l11 p11 (t ) = p11 = m (1 - p00 (t )) - mp11 (t ) = m - m ( p00 (t ) + p11 (t )) . Решения, естественно, совпадают с найденными. Листинг программы и полученные результаты при l = 4 m и = 7. l := 4m := P00 ( t) := l 7 e - ( l + m ) Чt l +m +m P01 ( t) := -l e - ( l + m ) Чt l +l +m На рис. 7.2 представлены графики этих функций. 124 7.4. Модель многоканальной СМО с отказами 0.8 P00 ( t ) 0.6 P01 ( t ) 0.4 0.2 0.2 0.4 0.6 0.8 t Рис. 7.2 7.4. Модель многоканальной СМО с отказами Математическая модель Пусть есть n — каналов обслуживания, и заявки поступают на обслуживание, если свободен хотя бы один канал. В против‑ ном случае заявка получает отказ. Система может находиться в n + 1 состоянии: 0 — все каналы свободны; 1 — свободен один канал; 2 — свободны два канала; … n — свободны n‑каналов. В этом случае матрица плотностей переходов имеет вид м-l, i = 0 п l 0i = н l, i = 1 , п 0, i > 2 о м im, j = i - 1 п li j = н - ( l + im ) , j = i , п l, j = i + 1 о мnm, j = n - 1 ln j = н , о -nm, j = n где λ — параметр потока заявок, а μ — параметр обслуживания. 125 7. Непрерывные стохастические модели (Q‑схемы) Таким образом, ж -l з зm з0 з з0 L=з з з з0 з з и0 l - (l + m ) m  l - ( l + 2m ) m  l - ( l + 3m ) m  ј ј ј ј ј  ј ј nm  - ( l + (n - 1) m ) ц ч ч ч ч ч ч. ч ч ч ч -nm чш  l Система уравнений Колмогорова P ў (t ) = P (t ) L имеет вид: м piў0 (t ) = - pi 0 (t ) l + pi1 (t ) m, i = 0,1, ј, n п н piў j (t ) = lpi j -1 (t ) - ( l + m ) pi j (t ) + ( j + 1) mpi j +1 (t ) , i = 0, 1, ј, n, j = 1, 2, ј, n - 1 п ў о pin (t ) = pi n -1 (t ) l - nmpi n (t ) , i = 0,1, ј, n, с начальными условиями pi j ( 0 ) = 0, i № j , pi i ( 0 ) = 1, i = 0, 1, , n. Из структуры матрицы видно, что система разбивается на n +1 блок уравнений, в каждый из которых входят неизвестные одной строки. В случае системы с двумя каналами обслуживания имеем три состояния: l 0 ц ж -l ч з L = з m - (l + m ) l ч з 0 2m -2m чш и и м piў0 (t ) = - pi 0 (t ) l + pi1 (t ) m п н piў1 (t ) = lpi 0 (t ) - ( l + m ) pi 1 (t ) + 2mpi 2 (t ) , п ў о pi 2 (t ) = pi 1 (t ) l - 2mpi 2 (t ) получаем систему уравнений Колмогорова, которая разбивает‑ ся на три блока: i = 0, 1, 3. 126 7.4. Модель многоканальной СМО с отказами Решение данной системы можно получить аналитически. Однако при большом количестве уравнений в системе реше‑ ние проще получить, используя уже известные по главе 2 про‑ цедуры пакета MathCAD. Ниже приведен листинг программы и результаты расчета при l = 4, m = 7 и pi i ( 0 ) = 1, pi j ( 0 ) = 0, i № j , см. рис. 7.3. -l Чp0 + m Чp1 й щ йж 1 ц щ к ъ D ( t , p) := к l Чp0 - ( l + m ) Чp1 + 2m Чp2 ъ P := rkfixed кз 0 ч , 0 , 1.1 , 100 , Dъ кз ч ъ к ъ ли 0 ш ы l Ч p1 - 2m Ч p2 л ы t := Pб0с P00 := Pб1с P01 := Pб2с P02 := Pб3с . P00 0.8 P01 0.6 P02 0.4 0.2 0.2 0.4 0.6 0.8 t Рис. 7.3 Стационарные вероятности получим из системы м- pi 0 l + pi1m = 0 п нlpi 0 - ( l + m ) pi 1 + 2mpi 2 = 0 п о pi 1l - 2mpi 2 = 0. с учетом, что pi 0 + pi 1 + pi 2 = 0. Из первого уравнения следует m pi 0 = pi1 . Подставляя в третье уравнение, получим l 127 7. Непрерывные стохастические модели (Q‑схемы) lpi 1 = 2m - 2m 2 2ml pi 1 - 2mpi 1 и pi 1 = 2 . l l + 2m 2 + 2ml Тогда pi 0 = pi 2 = 1 - 2m 2 , l 2 + 2m 2 + 2ml 2m 2 2ml l2 . = l 2 + 2m 2 + 2ml l 2 + 2m 2 + 2ml l 2 + 2m 2 + 2ml Из полученных результатов следует, что стационарные ве‑ роятности не зависят от i, т. е. от того, в каком состоянии нахо‑ дилась система изначально. Эта ситуация есть отражением того факта, что однородный и транзитивный марковский процесс является эргодическим [11], т. е. существуют пределы lim pi j (t0 ,t0 + Dt ) = lim pi j ( Dt ) = p j > 0. Эргодический марковский Dt ® 0 Dt ® 0 процесс с конечным числом состояний i = 1, 2, , n всегда об‑ ладает стационарным режимом. В этом случае в качестве харак‑ теристик процесса рассматривают вероятности pi (t ) нахожде‑ ния системы в i -ом состоянии и уравнения Колмогорова для определения этих функций имеют вид n n n n k =1 k №i k №i k №i piў (t ) = е pk (t ) l k i = е pk (t ) l k i + pi (t ) li i = е pk (t ) l k i - pi (t ) е li k c начальными условиями pi ( 0 ) = pi0 . Стационарные решения на‑ n n k №i k №i ходим из системы однородных уравнений е pk lk i - pi е li k = 0 n с учетом нормировки е pk = 1. k =1 Для многоканальной СМО с отказами эта система имеет вид м p0ў (t ) = - p0 (t ) l + p1 (t ) m п н piў (t ) = lpi -1 (t ) - ( l + im ) pi (t ) + (i + 1) mpi +1 (t ) , i = 1, 2,ј, n - 1 п ў о pn (t ) = pn -1 (t ) l - nmpn (t ) . 128 7.5. Модель многоканальной СМО с очередью 7.5. Модель многоканальной СМО с очередью Математическая модель Пусть система имеет n каналов обслуживания и в очереди может быть не более m заявок, т. е. заявка отклоняется, если в очереди уже m заявок. Максимальное число заявок в системе n + m , а число состояний системы n + m +1. Система уравнений Колмогорова для pi (t ) имеет вид ж p0ў (t ) = - p0 (t ) l + p1 (t ) m з з piў (t ) = lpi -1 (t ) - ( l + im ) pi (t ) + (i + 1) mpi +1 (t ) , i = 1, 2, ј, n - 1 з pў (t ) = lp (t ) - ( l + im ) p (t ) + nmp (t ) , i = n i -1 i i +1 з i з piў (t ) = lpi -1 (t ) - ( l + nm ) pi (t ) + nmpi +1 (t ) , i = n + 1, ј, n + m - 1 з з pў (t ) = p n + m -1 (t ) l - nmpn + m (t ) . и n+m Рассмотрим решение этой системы для двухканальной си‑ стемы с накопителем емкостью в 2 заявки. В этом случае систе‑ ма дифференциальных уравнений состоит из пяти уравнений ж p0ў (t ) = - p0 (t ) l + p1 (t ) m з з p1ў (t ) = lp0 (t ) - ( l + m ) p1 (t ) + 2mp2 (t ) з pў (t ) = lp (t ) - ( l + 2m ) p (t ) + 2mp (t ) 1 2 3 з 2 з p3ў (t ) = lp2 (t ) - ( l + 2m ) p3 (t ) + 2mp4 (t ) з з pў (t ) = p (t ) l - 2mp (t ) . 4 5 и 5 Ниже приведены листинг программы, результаты для на‑ чального состояния — одна заявка в системе обслуживания и графики решения (рис. 7.4). 129 7. Непрерывные стохастические модели (Q‑схемы) - l Ч p0 + m Ч p1 йк ъщ щ йж 0 ц к l Чp0 - ( l + m ) Чp1 + 2m Чp2 ъ з к ъ 1ч к ъ l Ч p1 - ( l + 2m ) Ч p2 + 2m Ч p3 ъ P := rkfixed кз 0 ч , 0 , 1.1 , 100 , Dъ D ( t , p) := к ъ кз 0 ч к l Чp2 - ( l + 2m ) Чp3 + 2m Чp4 ъ ъ кз ч к ъ з ч ъ к l Ч p l 2 m ( + ) Ч p 2 + m Ч p 4 5ъ к 3 ли 0 ш ы к ъ l Ч p4 - 2m Ч p5 л ы t := Pб0с P0 := Pб1с P1 := Pб2с P2 := Pб3с P3 := Pб4с P4 := Pб5с P5 := Pб6с P0 P1 P2 0.8 0.6 P3 P4 0.4 P5 0.2 0.5 1 1.5 t Рис. 7.4 Стационарные вероятности в общем случае находятся по формулам 1 p0 = i i n+m n n ж l ц n -1 1 ж l ц + е з ч е з ч i = n n! и nm ш i = 0 i! и nm ш i ж l ц pi = p0 з ч , 1 Ј i Ј n i! и nm ш i pi = p0 130 nn ж l ц з ч , i > n. n! и nm ш Приложение Встроенные функции MathCAD 1. К главам 2 и 7 Численное решение обыкновенных дифференциальных уравнений и систем. Процедура — функция Odesolve ([vector], x, b, [intvls]): vector — вектор названий функций (используется только для систем), х — имя переменной, b — конечная точка интервала интегрирования (должна быть больше точки, указанной в начальных условиях, intvls — целое число интервалов, используемых для интер‑ поляции полученного решения. Odesolve возвращает решение в виде вектор-функции неза‑ висимой переменной, значения которой в расчетных точках ис‑ пользуются для интерполяции функцией ispline. Odesolve предполагает, что старшие производные входят в систему линейно. Odesolve использует по умолчанию (Adams/BDF), для реше‑ ния нежестких систем выбирается метод Адамса — Бэшворда (Adams), а для решения жестких систем — метод формул диф‑ ференцирования назад (BDF) [13]. Однако если в блоке реше‑ ния кликнуть правой клавишей на названии Odesolve, то в по‑ явившемся меню, кроме Adams/BDF, можно выбрать еще три позиции: 131 Приложение 1) Fixed — соответствует использованию процедуры-функ‑ ции rkfixid; 2) Adaptive — соответствует использованию процедурыфункции rkadapt; 3) Radau — соответствует использованию процедуры-функ‑ ции Radou. Процедура-функция rkfixid (y, x1, x2, n, D): y — вектор начальных значений искомых функций, x1, x2 — начальное и конечное значение области изменения переменной, n — число точек, в которых ищется решение, D — вектор-функция правых частей cистемы уравнений. Процедура rkfixid использует пятиточечный метод Рунге — Кутта с фиксированным шагом. Процедура rkfixid возвращает матрицу, где первый столбец есть набор точек, в которых было найдено решение, а после‑ дующие столбцы представляют значения искомых функций в этих точках. Процедура rkfixid применима и в случае нелинейного вхож‑ дения старшей производной в уравнение. Процедура-функция rkadapt (y, x1, x2, acc, D, kmax, s): y — вектор начальных значений искомых функций, x1, x2 — начальное и конечное значение области изменения переменной, acc — точность интегрирования, D — вектор-функция правых частей cистемы уравнений, kmax — максимальное число точек, s — минимально допустимый интервал между точками. Процедура rkadapt использует метод Рунге — Кутта с пере‑ менным шагом. Процедура-функция Radou (init, x1, x2, intvls, D, [J], [M], [tol]) применяется для решения жестких систем с алгебраиче‑ скими ограничениями. 132 Встроенные функции MathCAD Примеры применения процедур Odesolve и rkfixid 1. Решить задачу Коши y ўў ( x ) + xy ў ( x ) = sin ( x ) , y ( 0 ) = 0, y ў ( 0 ) = 1. Given 0 y' ( 0) 0 y ( 0) y'' ( x) + xЧy' ( x) – sin( x) 1 y := odesolve( x , 10) 3 2 y( x) 1 2 4 6 8 10 x пмy1ў (t ) = ty1 (t ) - y2 (t ) + t 2. Решить систему н оп y2ў (t ) = 3 y1 (t ) + y2 (t ) , y1 ( 0 ) = 1, y2 ( 0 ) = 2. а) с использованием процедуры Odesolve: Given d y1 ( t) dt d y2 ( t) dt t Чy1 ( t) – y2 ( t) – t y1 ( 0) y1 ( t) + y2 ( t) y2 ( 0) йж y1 ц щ ж y1 ц з ч := Odesolve кз ч , t , 2, 5ъ ли y2 ш ы и y2 ш 1 4 y1 ( t ) y2 ( t ) 2 -2 -4 0.5 1 1.5 2 t 133 Приложение б) с использованием процедуры rkfixid: Given й t Чy0( t) – y1( t) – t щ ъ R := rkfixed йлжи 1 цш , 0, 2, 100,Dщы л y0( t) + y1( t) ы D ( t , y) := к y1 := Rб1с y2 := Rб2с 4 y1 y2 2 -2 -4 0.5 1 1.5 2 Rб0с 2. К главе 3 Численные методы нахождения безусловных и условных экс‑ тремумов. Процедуры minimize и maximize. Процедура-функция minimize (f, var1, var2, …). f(var1, var2, …) — минимизируемая функция. Задается вне блока решения. var1, var2, … — переменные минимизируемой функции. Перед использованием вне блока решения необходимо за‑ дать приближенные начальные значения var1, var2, … . Процедура-функция maximize (f, var1, var2, …). Максимизирует значение f(var1, var2, …). Обе процедуры возвращают значения var1, var2, … в найден‑ ных точках минимума или максимума. 134 Встроенные функции MathCAD Ограничения на область, в которой ищется минимум или максимум, задаются в блоке решения. Процедуры используют для решения различные методы из целой группы методов в за‑ висимости от того, является ли задача линейной или нет, а так‑ же от характера ограничений. Рассмотрим примеры. а) Найти минимум функции z = x 2 + 2 x + y 2 . Листинг программы. z ( x , y) := x2 + 2x + y 2 x := 0.5 y := 0.5 Given ж r := minimize ( z , x , y) r = з и -1 ц. ч -1.686 ґ 10- 9 ш Найдены координаты точки минимума (–1, 0). б) Найти минимум функции z = x 2 + 2 x + y 2 при условии y = x . Листинг программы. z ( x , y) := x2 + 2x + y2 x := 0.5 y := 0.5 Given x -0.5 ц . y r := Minimize ( z , x , y) r = ж и -0.5 ш 3. К главе 6 Процедуры обработки экспериментальных данных. 1. linfit (vx, vy, F) производит подгонку с использованием функции y ( x ) = k1 f1 ( x ) + k2 f 2 ( x ) +  + kn fn ( x ) : vx — массив значений переменной; vy — массив значений функции y ( x ); ж f1 ( x ) ц ч з f2 ( x ) ч . F — вектор-функция зз  ч чч зз и fn ( x ) ш 135 Приложение Процедура linfit (vx, vy, F) возвращает вектор коэффици‑ ж k1 ц з ч k2 ентов з ч . з  ч зз чч и kn ш 2. genfit (vx, vy, vs, F) производит подгонку с использовани‑ ем функции f ( x, k ) , где k — вектор параметров: vs — начальные значения параметров; F — вектор, содержащий функцию f ( x, k ) и ее производные по параметрам. Процедура возвращает значения вектора параметров. 3. regress (vx, vy, n) — полиномиальная регрессия порядка n: n — степень полинома. Возвращает вектор коэффициентов vs, запрашиваемый про‑ цедурой interp (vx, vy, vs, x). Примеры использования этих процедур. Пусть VX = ( -1, 0, 1, 2, 3, 4 ) ,VY = (1.1, 0, 0.9, 3.5, 7, 12 ) . 1. Простая и квадратичная регрессии с использованием про‑ цедуры linfit (vx, vy, F). Листинг программы (простая регрессия). T T ж 1ц vx := ( 1 – 1 0 2 3 4) vy := ( 1 0 0.8 3 7 13) F1 ( x) := з ч и xш 0.465 ц r := linfit ( vx , vy , F1) = ж и 2.446 ш Квадратичная регрессия. ж 1ц ж -0.026 ц з ч x r := linfit ( vx , vy , F2) = з 0.239 ч F2 ( x) := з 2ч ч з и 0.736 ш иx ш 136 Встроенные функции MathCAD 2. Квадратичная регрессия с использованием процедуры genfit (vx, vy, vs, F). Листинг программы. ж b 0 + b 1 Чx + b 2 x2 ц з ч 1 T з ч vs := ( 0 1 1) F3 ( x , b ) := з ч x з ч x2 и ш ж -0.026 ц r := genfit ( vx , vy , vs , F3) = з 0.239 ч з ч и 0.736 ш Результаты полностью совпадают с полученными выше. 3. Подгонка с помощью функции f ( x, k ) = k0e - x + k0 k1 x и процедуры genfit (vx, vy, vs, F). Листинг программы. 1 -x vs := ж ц f ( x , k ) := k 0 Чe + k 0 Чk 1 Чx и 1ш ж k 0 Ч e- x + k 0 Ч k 1 Ч x ц з ч з ч r := genfit ( vx , vy , vs , F4) = ж 0.824 ц -x F4 ( x , k ) := e + k Ч x 1 и 3.19 ш з ч з ч k 0 Чx и ш Функции математической статистики 1. mean (X) — выборочное среднее. Вычисляется по форму‑ 1 n -1 ле е xi . n i =0 2. Var (X) — исправленная выборочная дисперсия. Вычис‑ 2 1 n -1 ляется по формуле е ( xi - mean ( X )) . n - 1 i =0 2 1 n -1 3. var ( X ) = е ( xi - mean ( X )) — выборочная дисперсия. n i =0 137 Приложение 4. Stdev ( X ) = Var ( X ). 5. stdev ( X ) = var ( X ). 1 n -1 6. cvar ( X ,Y ) = е ( xi - mean ( X )) ( yi - mean (Y )) — выбороч‑ n i =0 ная ковариация. c var( X ,Y ) 7. corr( x, y ) = — коэффициент корреляции. stdev( X )stdev(Y ) 8. int ercept(X,Y) — коэффициент b0 в модели регрессии y = b0 + b1 x + e. 9. slope(X,Y) — коэффициент b1 в модели регрессии y = b0 + b1 x + e. 10. stderr(X,Y) = 2 1 n -1 yi - int ercept ( X ,Y ) - slope( X ,Y )xi ) — ( е n - 2 i =0 оценка среднеквадратичной ошибки σ. 11. qt ( p, d ) — квантиль t p (d ) распределения Стьюдента с d степенями свободы. 12. qF ( p, d1, d 2 ) — квантиль F p (d1, d 2 ) распределения Фише‑ ра с d1 и d 2 степенями свободы. 13. qchisq ( p, d ) — квантиль c2p (d ) распределения хи-квадрат с d степенями свободы. 138 Библиографический список 1. 2. 3. 4. 5. 6. 7. 8. Тихонов А. Н. Уравнения математической физики : учеб‑ ное пособие для университетов / А. Н. Тихонов, А. А. Са‑ марский. — Москва : Наука, Гл. ред. физ.-мат. лит., 1966. — 718 с. Советов Б. Я. Моделирование систем : учебн. для ву‑ зов / Б. Я. Советов, С. А. Яковлев. — Москва : Высш. шк., 2001. — 343 с. Самарский А. А. Математическое моделирование. Идеи. Методы. Примеры : моногр. / А. А. Самарский, А. П. Ми‑ хайлов. — Москва : ФИЗМАТЛИТ, 2002. — 320 с. Математическое моделирование : пер. с англ. / под ред. Дж. Эндрюс, Р. Мак-Лоун. — Москва : Мир, 1979. — 276 с. Моисеев Н. Н. Методы оптимизации: моногр. / Н. Н. Мо‑ исеев, Ю. П. Иванилов, Е. М. Столярова. — Москва : На‑ ука, Гл. ред. физ.-мат. лит., 1978. — 352 с. Схрейвер А. Теория линейного и целочисленного про‑ граммирования: моногр. В 2‑х т. Т. 2 / А. Схрейвер. — Москва : Мир, 1991. — 342 с. Лоскутов А. Ю. Введение в синэргетику : учебн. руковод‑ ство / А. Ю. Лоскутов, А. С. Михайлов. — Москва : Нау‑ ка, Гл. ред. физ.-мат. лит., 1990. — 272 с. Тоффоли Т. Машины клеточных автоматов / Т. Тоффо‑ ли, Н. Марголус ; пер. с англ. под ред. Б. М. Баталова. — Москва : Мир, 1991. — 280 с. 139 Библиографический список 9. 10. 11. 12. 13. 140 Себер Дж. Линейный регрессионный анализ / Дж. Се‑ бер ; пер. с англ. под ред. М. Б. Малютова. — Москва : Мир, 1980. — 456 с. Дрейпер Н. Прикладной регрессионный анализ / Н. Дрейпер, Г. Смит ; пер. с англ. под ред. М. Саит-Ах‑ метова. — Москва : Изд. дом «Вильямс», 2007. — 912 с. Боровков А. А. Теория вероятностей : учебное пособие для вузов / А. А. Боровков. — Москва : Наука, Гл. ред. физ.-мат. лит., 1986. — 432 с. Вентцель Е. С. Теория случайных процессов и ее инже‑ нерные приложения : моногр. / Е. С. Вентцель, Л. А. Ов‑ чаров. — Москва : Наука, Гл. ред. физ.-мат. лит., 1991. — 384 с. Современные численные методы решения обыкновен‑ ных дифференциальных уравнений / под ред. Дж. Хол‑ ла и Дж. Уатта ; пер. с англ. под ред. А. Д. Горбунова. — Москва : Мир, 1979. — 312 с. Оглавление Предисловие......................................................................3 1. Понятие математической модели. Некоторые типовые схемы математического моделирования......5 1.1. Формальная математическая модель объекта моделирования.......................................................5 1.2. Общая схема создания математической модели и работы с ней...........................................6 1.3. Типовые математические схемы моделирования.......................................................7 2. Непрерывно-детерминированные схемы математического моделирования (D‑схемы)..............8 2.1. Эволюционные модели.........................................8 2.2. Некоторые модели математической физики......13 2.3. Задачи...................................................................22 3. Модели оптимизации.................................................23 3.1. Условный и безусловный экстремум дифференцируемой функции..............................23 3.2. Задачи линейного программирования................28 3.3. Нелинейное программирование..........................45 3.4. Выпуклое программирование..............................52 3.4. Задачи...................................................................53 141 Оглавление 4. Дискретно-детерминированные схемы (F‑схемы)....55 4.1. Конечные автоматы.............................................55 4.2. Клеточные автоматы............................................60 4.3. Задачи...................................................................65 5. Дисктетно-стохастические схемы (Р‑схемы)............67 5.1. Вероятностные автоматы.....................................67 5.2. Задачи...................................................................74 6. Регрессионные модели...............................................75 6.1. Модель одномерной линейной регрессии..........75 6.2. Модель множественной регрессии......................90 6.3. Одномерные полиномиальные модели............. 102 6.4. Планирование регрессионного эксперимента.... 109 6.5. Задачи................................................................. 114 7. Непрерывные стохастические модели (Q‑схемы)... 117 7.1. Системы массового обслуживания.................... 117 7.2. Марковские процессы....................................... 119 7.3. Модель одноканальной СМО с отказами......... 123 7.4. Модель многоканальной СМО с отказами....... 125 7.5. Модель многоканальной СМО с очередью....... 129 Приложение................................................................... 131 Библиографический список......................................... 139 142 Учебное издание Мохрачева Людмила Павловна Типовые математические схемы моделирования. Примеры и задачи Редактор О. В. Климова Корректор А. А. Трофимова Верстка О. П. Игнатьевой Подписано в печать 17.04.2018. Формат 60×84/16. Бумага офсетная. Цифровая печать. Усл. печ. л. 8,4. Уч.-изд. л. 4,9. Тираж 50 экз. Заказ 102 Издательство Уральского университета Редакционно-издательский отдел ИПЦ УрФУ 620049, Екатеринбург, ул. С. Ковалевской, 5 Тел.: +7 (343) 375-48-25, 375-46-85, 374-19-41 E-mail: [email protected] Отпечатано в Издательско-полиграфическом центре УрФУ 620083, Екатеринбург, ул. Тургенева, 4 Тел.: +7 (343) 358-93-06, 350-58-20, 350-90-13 Факс: +7 (343) 358-93-06 http://print.urfu.ru 9 785799 623623