Аннотация В первой части работы говорится об особенностях моделирования различного рода динамических процессов и приводится пример работы с одной из таких моделей - Кэйбабского плато: от постановки задачи до анализа ее результатов. Во второй части подробно описывается версия языка DYNAMO, являющаяся входным языком компилятора для персональных ЭВМ. Компилятор DYNAMO реализован как препроцессор над C (используется версия Quick C 2.0 фирмы Microsoft). Часть I. Моделирование динамических процессов. Введение Моделирование в широком понимании - это естественный и едва ли не единственный метод познания природы. И действительно, все наши знания об окружающем нас мире представляются в виде самых разнообразных моделей: интуитивных, словесных, механических, электрических, математических. Различать и классифицировать модели человек начал сравнительно недавно, но, несомненно, бессознательное моделирование - ровесник разумного человека. Возьмем пример сегодняшнего дня. Особенно остро стоит сейчас вопрос о создании комплексных моделей, описывающих экономические, географические, метеорологические, гидрогеологические, демографические и многие другие процессы, по возможности учитывающие многообразные сложные связи между ними. Без таких моделей, видимо, не представляется возможной реализация долгосрочных социально-экономических программ, проведение крупномасштабных мероприятий в области охраны окружающей среды, освоение новых хозяйственных районов и составление глобальных планов развития экономики. Все эти проблемы не новы, но в наше время комплексное решение их нельзя назвать просто желательным - это насущная потребность сегодняшнего дня, это вопрос о том, как будем выглядеть мы сами, наши города, наши заводы, наша природа, короче говоря, каким будет завтрашний день нашей планеты. Эта естественная заинтересованность общества находит сейчас отражение во многих государственных актах. Для успешного претворения этих проектов в жизнь нужна прочная научная основа. Естественно, что модели, которые создает человек, совершенствуются вместе с методами исследований, с арсеналом средств для моделирования, которыми располагает наука. Всякая отрасль науки неизбежно переживает в своем развитии три периода - период накопления сведений, период обобщения и установления эмпирических закономерностей, на базе которых создаются модели на уровне причинных связей, и, наконец, период, который мы назовем периодом математизации. Этому периоду свойственно установление четких количественных зависимостей и создание именно математических моделей. Наступление этого периода закономерно и необходимо для дальнейшего развития науки. Нельзя забывать, что постановка натурного, физического эксперимента с современными сложными системами иногда бывает чрезвычайно затруднительна, а порой и невозможна по многим причинам. Так, причина и следствие часто бывают разнесены во времени и пространстве. Особенно часто это наблюдается, например, в явлениях, происходящих в биосфере Земли и связанных с внешней средой обитания человека. Известно, что вредное действие некоторых инсектицидов (ДДТ и тому подобные), которые довольно интенсивно применялись в сельском хозяйстве во всем мире, может проявится только через десятки лет. Примером разнесения причины и следствия на расстояние может служить факт накопления свинца во льдах Гренландии как результат увеличения использования этого материала в Европе и сбрасывания его в море в виде промышленных отходов. Нельзя не отметить и то, что первая непосредственная реакция на внешние воздействия в сложных системах зачастую обратна тем изменениям, которые произойдут через более продолжительное время. Например, некоторые сыворотки, применяемые для противоэпидемических прививок, вызывают вначале ухудшение самочувствия, хотя окончательных итог их действия - повышения сопротивляемости организма к эпидемическим инфекциям. Использование математических моделей позволяет избегать дорогостоящих натурных экспериментов, сокращает время на проверку гипотез, открывает возможность многократного повторения опытов и устраняет многие другие трудности. Конечно, и математические модели не являются панацеей от всех бед, но именно они позволяют максимально использовать имеющиеся сведения и создавать разумные гипотезы, т.е. являются весьма действенным инструментом познания окружающего мира, средством прогнозирования и надежного управления. Что такое модель и имитационная система Под моделью мы будем понимать относительную истину, отражающую определенные особенности изучаемых явлений. Модель отражает наши представления о реальном мире. Вместе с изменением наших представлений о реальном мире по мере развития наших знаний должны меняться и описывающие его модели. Наличие гармоничной системы моделей, где более общие, но сложные модели органично включают в себя более простые в качестве частных случаев, отражает глубину наших представлений о сущности изучаемых явлений. Во многих областях науки такие системы моделей уже существуют, в других областях еще предстоит создать и развить их. От всех моделей математические модели отличаются тем, что средством описания моделей и изучения их поведения является формально-логический аппарат математики. Отсюда следует важнейшее преимущество - широкая возможность количественного анализа моделей с помощью современных математических методов. Другим важным преимуществом математических моделей является универсальность языка математики, возможность использовать одни и те же модели для исследования физически различных систем. Так, уравнения движения материальной точки в поле тяготения представляют собой математическую модель чрезвычайно широкого класса реальных явлений. Эта модель описывает как движение планет Солнечной системы, так и полет ракеты в космическом пространстве. Заметим при этом, что адекватного физического объекта для понятия "материальная точка", т.е. объекта, обладающего конечной массой, сосредоточенной в нулевом объеме, не существует. Еще одним полезным свойством математических моделей является возможность получать результаты, относящиеся не к отдельной конкретной реализации, соответствующей определенным начальным данным и фиксированным значениям параметров исследуемой системы, а сразу для целого множества возможных поведений системы. Это обстоятельство облегчает качественный анализ математической модели, что проявляется особенно ярко в тех случаях, когда удается получить описание поведения системы с помощью аналитических выражений, формул. Модель или совокупность моделей, описывающих данное явление, вместе со средствами анализа их поведения образуют имитационную систему. Примером имитационной системы, позволяющей определять основные аэродинамические характеристики реального самолета, может служить комплекс, состоящий из макета самолета, аэродинамической трубы и аппаратуры для измерения и обработки результатов. К сожалению, для большинства математических моделей, описывающих современные сложные системы получение аналитических результатов невозможно. Поэтому приходится прибегать к использованию алгоритмов численного анализа, что, в свою очередь, делает необходимым применение в имитационных системах электронных вычислительных машин (ЭВМ). Таким образом, в настоящее время можно говорить о создании математических имитационных систем, использующих ЭВМ. Но прежде чем приступить к рассмотрению различных аспектов использования ЭВМ в имитационных системах, рассмотрим методологию и методику построения математических моделей реальных явлений. Как строить математические модели. Этап выработки концепций. Сейчас уже накоплен определенный опыт построения математических моделей в таких нетрадиционных для приложения математики областях, как экономика, биология, медицина, организация промышленного производства и т.п. Всякое моделирование осуществляется в два этапа. Первый этап состоит в изучении всей комплексной проблемы, для решения которой (в полном объеме или в какой-то ее части) мы намереваемся применить имитационную систему. Назовем условно этот этап этапом выработки концепций. Второй этап - этап реализации выработанных на первом этапе концепций, которые представляются в виде математических соотношений между параметрами системы (рис. 1). Реальный ===> Восприятие реального ======> Модель - реализация мир мира и сведения о нем идей разработчиков I у разработчика модели II модели этап этап "концептуализация" "реализация" Рис. 1. Два этапа процесса моделирования. Этап выработки концепций обычно включает в себя следующие шаги: а) внутри изучаемой проблемы выделяется задача, которая кажется подходящей для исследования методами математического моделирования. Следует отметить, что отнюдь не все задачи являются таковыми. Причинами этого может быть, например, недостаточная изученность отдельных частей исследуемого явления. Результатом первого шага должно явиться словесное описание выбранной задачи; б) формулируются те специфические вопросы, на которые должен быть дан ответ после моделирования; в) выбирается длина интервала времени, на котором будет проводиться моделирование. Выбранный интервал времени должен быть достаточным для того, чтобы успели проявиться все характерные черты данного явления. Выбор неоправданно большого интервала времени чаще всего затрудняет анализ поведения системы, особенно в тех случаях, когда реализация модели связана с применением ЭВМ; г) формулируются ограничения на рассматриваемую систему и принимается решение о степени детализации модели, т.е. определяется количество параметров, включаемых в модель, число и сложность связей между этими параметрами, виды управляющих воздействий на систему. Такое решение невозможно без решения вопросов, поставленных в пунктах "б" и "в". Действительно, те специфические вопросы, на которые мы хотим получить ответ с помощью математического моделирования, будут определять и степень детализации модели. Выбор временного интервала точно так же существенно влияет на возможность реализации модели. С одной стороны, выбор короткого интервала времени моделирования системы часто позволяет отбрасывать связи, которые не оказывают заметного действия на поведение системы на рассматриваемом интервале времени. С другой стороны, увеличение временного интервала может привести к включению в модель новых факторов, действие которых со временем оказывается существенным для модели. Это, в свою очередь, приводит к появлению в модели новых параметров и новых, быть может, сложных зависимостей. В конце концов модель окажется очень громоздкой, что затруднит и даже сделает невозможным реализацию модели на ЭВМ; д) выбирается уровень агрегирования системы. Под агрегированием мы понимаем замену нескольких параметров системы одним параметром. Понятие уровня агрегирования тесно связано с понятием степени детализации, однако не тождественно ему. Если степень детализации в основном зависит от количества факторов как внутренних, так и внешних, включенных в рассмотрение, то уровень агрегирования определяется возможностью объединять в один фактор группы различных факторов, сходным образом действующих на систему. Это, естественно, упрощает модель. Так как простые модели, в общем, предпочтительнее, следует выбирать максимальный уровень агрегирования, при котором еще возможно отвечать на вопросы, сформулированные в пункте "б"; е) составляется диаграмма связей зависимостей, т.е. ориентированный граф, в вершинах которого указаны названия параметров модели, а направление дуги, например, от вершины А к вершине В означает, что изменение параметра А вызывает изменение параметра В. Иногда бывает полезным помечать дугу значком "+" или "-" в зависимости от того, приводит ли увеличение параметра А к увеличению или к уменьшению В. Таким образом, получается легко обозримое изображение структуры моделируемой системы. Не следует рассматривать первый этап моделирования лишь как подготовительную работу для перехода к следующему этапу - реализации модели на ЭВМ. Выполнение перечисленных выше шагов прежде всего ведет к четкому пониманию поставленной задачи, выявляет пробелы в наших знаниях о системе, заставляет обсуждать и принимать те или иные гипотезы. Уже на этом этапе широко используются эксперименты, для проведения которых привлекаются специалисты в тех конкретных областях знания, которые относятся к модели. Вряд ли будет успешным моделирование, если к нему уже на первом этапе работы не будет привлечен будущий пользователь - "заказчик" модели. Использование модели начинается задолго до того, как будет завершена работа по созданию имитационной системы! Как правило, анализ диаграммы причинных зависимостей позволяет обнаружить существенные качественные эффекты в поведении системы. Последнее обстоятельство часто приводит к необходимости повторения всех или части шагов этапа выработки концепций и построению новой диаграммы причинных связей. Таким образом, первый этап носит итерационный характер, и изложенную нами последовательность действий на этом этапе следует рассматривать как идеальную. Проиллюстрируем первый этап построения модели следующим при- мером. Модель развития популяции оленей на Кэйбабском плато. Проблема строгого количественного изучения жизни сообществ животных, их взаимодействие с окружающей средой и сообществами других животных является весьма актуальной, поскольку без детального изучения этих вопросов нельзя предпринимать эффективных экономических мер по охране природы. В качестве подходящей для математического моделирования задачи была рассмотрена модель динамики численности оленей и хищников на Кэйбабском плато. Особенностью этого плато является его относительная изолированность от остального мира. Это обстоятельство позволяет ограничиться сравнительно небольшим числом внешних факторов, воздействующих на систему. Важным моментом для моделирования является наличие достаточно достоверных сведений о поголовье оленей и хищников на Кэйбабском плато в течении длительного периода времени (статистика велась с 1905 г.). Первый этап моделирования начнем со словесного описания системы: число оленей растет экспоненциально при отсутствии хищников. Увеличение числа оленей ведет к уменьшению и разрушению кормовой базы - истощаются и не успевают восстановиться в полной мере пастбища. В результате начинается вымирание оленей из-за недостатка пищи и численность популяции оленей уменьшается, приближаясь к некоторому равновесному уровню. Поставим пока скромную задачу получить характерные кривые изменения численности оленей на плато в зависимости от количества пищи на нем и числа хищников. Реальная кривая изменения численности оленей на плато приведена на рис. 2. Рис. 2. Изменение популяции оленей на Кэйбабском плато по годам. Временной интервал моделирования - приблизительно 50 лет. Ограничения на систему и степень детализации модели определяются условиями поставленной задачи. Основным ограничением системы является один и тот же постоянный район обитания (ареал) оленей. Для решения поставленной задачи (по крайнем мере в первом приближении) представляется возможным ограничиться в системе рассмотрением двух основных взаимодействий: связи численности волков и оленей (отношения хищник - жертва) и зависимости популяции оленей от кормовой базы. Ограниченность поставленной задачи и длительный временной интервал позволяют нам выбрать при моделировании достаточно высокий уровень агрегирования. Мы не будем различать оленей по полу, возрасту и другим признакам, а будем описывать популяцию одним параметром - числом животных. Аналогичным образом характеризуется и популяция хищников. В модели предполагается, что корм равномерно распределен по всей площади и не делается различия корма по виду, качеству и количеству. Теперь можно приступить к построению диаграммы причинных связей (рис. 3). Основными параметрами модели являются число оленей - DP, число хищников - PP, количество пищи - F, площадь плато - AREA. Рис. 3. Диаграмма причинных связей для модели Кэйбабского плато. Естественно предположить, что количество оленей, уничтожаемых хищниками за единицу времени (DPR), определяется количеством хищником (PP) и числом оленей, которых может убить один хищник (DKR); последнее, в свою очередь, зависит от плотности оленей DD. Можно было бы пояснить смысл и других параметров, изображенных на остальной части диаграммы, и характер взаимоотношений между ними. Однако делать это нам кажется нецелесообразным, поскольку назначение диаграммы причинных связей как раз и состоит в наглядном изображении этих взаимоотношений. Читатель без труда сможет сам разобраться в предположениях относительно динамики количества пищи и ее влияния на изменение поголовья оленей через естественный прирост стада. При построении приведенной здесь диаграммы ее авторы учитывали мнения экспертов - специалистов. Таким образом, диаграмма может рассматриваться как некоторый аккумулятор знаний о флоре и фауне изучаемого района. Оставим пока рассмотрение данного примера и перейдем к изложению методологии второго этапа моделирования. Реализация и анализ модели Итог первого этапа моделирования представлен диаграммой причинных связей, которая, несмотря на свою содержательность, не является на самом деле строгим, количественным описанием моделируемого явления. В действительности для того, чтобы получить результаты исследований в виде числовых характеристик поведения системы и получить их логически непротиворечивым путем, необходимо применить математические методы анализа. Собственно диаграмма причинных связей не поддается непосредственному исследованию каким-либо математическим методом, но она служит удобным исходным материалом для формализованного описания модели на языке той или иной математической дисциплины. Процесс развития математики всегда был тесно связан с процессом познания и описания природы. С этой точки зрения появление дифференциального и интегрального исчислений, теории дифференциальных уравнений - все это можно рассматривать как разработку языка для описания и анализа моделей ньютоновской механики и динамики сплошных сред. Развитие новых физических моделей - появление общей и специальной теорий относительности, квантовой механики, статистической физики и др. - повлекло за собой соответствующее развитие математического аппарата: возникли функциональный анализ и теория вероятностей, модифицировался старый аппарат, например, в теории дифференциальных уравнений стали рассматриваться обобщенные решения. Трудности получения аналитических решений во многих задачах математической физики и инженерного дела привели к созданию дискретной математики и конечноразностных моделей, аппроксимирующих непрерывные. Дискретные методы анализа оказались не только эффективными, но и единственно возможными при моделировании явлений, носящих дискретный характер. В качестве примеров можно привести экономические задачи, дискретность которых определяется заранее установленными периодами планирования и отчетности. Во многих конкретных задачах приходится иметь дело с системами смешанного типа, где одни процессы, непрерывные во времени, взаимодействуют с другими процессами - дискретными во времени. Примером такого процесса может служить полет космической ракеты, в ходе которого коррекция траектории осуществляется включением и выключением двигателя. Зачастую многие непрерывные процессы без заметного ущерба для точности будущих результатов представляют в виде конечноразностных моделей. Такое представление дает возможность просчитывать эти модели непосредственно на ЭВМ. Выбор математического языка для реализации модели является одним из ключевых моментов в процессе моделирования. Выбор этот производится с одновременным учетом двух факторов: природы моделируемого явления и возможностей анализа модели (аналитического, машинного или человеко-машинного). Ранее уже говорилось о невозможности серьезного аналитического исследования моделей современных сложных систем. Поэтому при выборе языка описания таких систем необходимо заранее ориентироваться на использование ЭВМ. В связи с этим большое значение приобретает разработка языков моделирования. Сейчас известно много языков, специально предназначенных для простого и удобного описания большого числа параллельно протекающих во времени взаимодействующих процессов. Во всех языках моделирования предусмотрена возможность ввода и работы с информацией в виде графиков, таблиц, т.е. в том виде, в котором обычно представляются результаты экспертиз и экспериментов. Чтобы моделирование стало действительно эффективным методом познания, необходимо привлекать специалистов в той области, в которой ведется моделирование, не только в качестве экспертов, а непосредственно в качестве авторов моделей. Это обстоятельство накладывает естественное требование на языки моделирования - они должны быть просты в употреблении и легкодоступны для быстрого изучения специалистами - нематематиками. Наконец, определенные требования предъявляются к математическому обеспечению той вычислительной техники, которая применяется в имитационных системах. Оно должно содержать в виде стандартных подпрограмм и процедур решение большинства типовых математических задач, встречающихся при обсчете моделей в данной конкретной области. Например, для имитирования экономических процессов необходимо, видимо, иметь программы, реализующие расчеты межотраслевого баланса, решение задач линейного программирования, массового обслуживания и др. Для моделирования биологических систем, очевидно, полезны программы решения задач статистического, корреляционного и регрессионного анализа. Для имитации технических систем требуются программы расчета изделий на прочность, устойчивость, программы расчета теплового баланса и других характеристик. Кроме того, математическое обеспечение должно содержать большое количество сервисных программ, позволяющих пользователю имитационной системы легко вносить изменения по ходу моделирования и получать результаты в привычном для человека виде - в виде таблиц, графиков, сводок. Это необходимо не только для экономии рабочего времени - пользователь не должен чувствовать себя стесненным системой. Он должен видеть в ней помощника, с которым легко и приятно работать. Этот психологический аспект ни в коем случае нельзя забывать. Конечно, такие мощные универсальные алгоритмические языки, как PASCAL, MODULA, C почти в полной мере удовлетворяют требованиям, предъявляемым к языкам моделирования. Но эти языки весьма сложны для изучения неспециалистом. Достаточно сложным является и такой универсальный язык моделирования, как SIMULA. Проблема сочетания простоты и достаточной выразительности языка моделирования может быть решена путем создания серии проблемно-ориентированных языков, т.е. языков, приспособленных для описания моделей реальных систем в конкретной области. В следующем параграфе мы кратко опишем язык моделирования DYNAMO, который, по нашему мнению, является удачным компромиссом, удовлетворяющим противоречивым требованиям простоты и универсальности. Знакомство с этим языком понадобится читателю для понимания формального описания модели, которые мы рассмотрим в этой книге. Речь идет о модели динамики популяции оленей на Кэйбабском плато. Вводные замечания к описанию языка моделирования DYNAMO Алгоритмический язык DYNAMO получил в свое время большое распространение в США. Этот язык подробно изложен во второй части данной работы. Он возник из потребностей алгоритмизации имитационных математических моделей, описывающих функционирование производственно-сбытовых систем крупных фирм. Указанные модели разрабатывались в США Форрестером [3]. До недавнего времени эти модели не представляли большого интереса для российского читателя, поскольку функционирование производственно-сбытовых подразделений в бывшем Союзе определялось совсем другими внутренними механизмами, чем в США. Однако, как это часто бывает, алгоритмический язык, возникший из потребностей решения некоторого частного класса задач, оказывается в дальнейшем удобным средством для описания более широкого круга задач и процессов. Такая ситуация имела место и с DYNAMO. Терминология, употребляемая в DYNAMO, и его структура дает возможность пользователям в достаточно естественной форме представлять течение процесса, если этот процесс сам по себе хорошо описывается системой обыкновенных дифференциальных уравнений. Язык DYNAMO предназначается для специалистов, не имеющих глубокого математического образования, однако детально представляющих изучаемый процесс. Изучив DYNAMO, такие специалисты получают возможность воспроизводить на ЭВМ и анализировать интересующий их процесс. Для изучения DYNAMO нет необходимости, чтобы специалист был знаком с теорией обыкновенных дифференциальных уравнений, хотя язык DYNAMO с формальной точки зрения описывает именно процесс численного интегрирования системы обыкновенных дифференциальных уравнений. Это обстоятельство является достоинством DYNAMO, и именно оно определяет круг его пользователей. Конечно, с точки зрения специалиста по вычислительной математике, знакомого с универсальными алгоритмическими языками, DYNAMO неудобен и обладает весьма ограниченными средствами и возможностями. Тем не менее, как показывает практика, круг пользователей DYNAMO в достаточной мере широк. В США большое число научно-исследовательских групп использовало DYNAMO в качестве основного рабочего инструмента. Благодаря DYNAMO непосредственный доступ к ЭВМ получило большое число специалистов, весьма поверхностно знакомых с математикой, таких, как управленческий персонал фирм, представители биологических специальностей и т.д. Все это позволяет предполагать, что может вызвать определенный интерес у соответствующих специалистов и в нашей стране. В языке DYNAMO используются для описания параметров моделируемой системы четыре вида переменных: уровни, темпы, вспомогательные переменные, дополнительные переменные. Обычно в качестве уровней в модели фигурируют те факторы (параметры) системы, которые численно описывают состояние основных процессов в моделируемой системе и динамику которых мы хотим получить в результате. В дальнейшем мы увидим, например, что в качестве уровней в модели Кэйбабского плато рассмотрены два параметра: численность оленей в данный момент времени и количество пищи. Скорости изменения уровней обычно задаются в DYNAMO с помощью темпов - переменных второго вида, указывающих, на сколько единиц и в какую сторону изменёется за единицу времени данный уровень, если остальные факторы (от которых зависит изменение значение этого уровня) остаются неизменными. Отметим, что темпы могут зависеть как от значения уровня, изменение которого они определяют, так и от значений других уровней, других темпов, а также от вспомогательных переменных. Примером темпа может служить в модели Кэйбабского плато естественный прирост оленьего стада за единицу времени. Вспомогательные переменные - это обычно такие параметры системы, получение динамики изменения которых не является целью моделирования. Они, однако, отражают и учитывают влияние различных факторов, в том числе и внешних, по отношению к системе, на темпы, управляющие изменением уровней модели. Если разработчик модели может сразу написать закон зависимости некоторого темпа от определяющих его параметров, то тем не менее вспомогательные переменные могут появиться и в этом случае, поскольку на форму описания зависимости между параметрами модели в DYNAMO наложены определенные ограничения. Наконец, дополнительные переменные употребляются при выводе результатов моделирования в том случае, если интересующая авторов модели величина не встречается в модели как уровень, темп или вспомогательная переменная. Например, можно предположить, что в модели Кэйбабского плато потребуется выводить на печать величину, означающую среднее количество оленей, приходящихся на одного хищника. Поскольку эта величина не является параметром модели, для ее представления необходимо ввести в модель соответствующую дополнительную переменную. Кроме переменных указанных четырех видов, во всякой DYNAMO-программе присутствует особая переменная - время. Это дискретная переменная, выбор единицы измерения которой, так же как и выбор временного шага (обозначаемого в DYNAMO символом DT), осуществляется разработчиком модели. в DYNAMO принято обозначать текущий момент времени буквой К, предшествующий - буквой J, будущий - L. Кроме переменных, в DYNAMO могут использоваться числа, константы. Основными символами языка DYNAMO являются буквы (латинского алфавита), цифры и ограничители. Буквы используются для образования идентификаторов - названий переменных и констант системы, из цифр составляются числа, а ограничители - это знаки операций: +, -, / (разделить), левая и правая скобки, = (знак равенства), десятичная точка, запятая, звездочка (она используется при описании переменных с индексами, которые могут составлять, например, содержание одновходовой таблицы или некоторого другого массива информации). Программа на языке DYNAMO представляет собой совокупность уравнений, описывающих взаимную зависимость параметров моделируемого процесса. Процесс выполнения программы заключается в вычислении по значениям величин, характеризующих динамический процесс в предыдущий момент времени, новых значений этих величин в последующий момент времени. Предыдущие моменты времени, как уже указывалось, обозначаются в DYNAMO символом J (или JK, если соответствующая величина является темпом). Соответствующие величины могут встречаться только в правых частях уравнений. Последующие моменты времени обозначаются в DYNAMO символом К (или KL для темпов). В заголовке программы указывается шифр программы, величина шага по времени, значение конца интервала времени. Вычисления величин, фигурирующих в левых частях уравнений программы, проводятся последовательно для всех моментов времени, начиная с нулевого, с шагом DT до момента времени, соответствующего концу интервала. На каждом временном шаге вычисления производятся в следующем порядке: сначала вычисляются все уровни в левых частях уравнений, относящиеся к моменту К, затем вычисляются вспомогательные переменные для момента К, после чего вычисляются все темпы на интервале KL. Если на данном временном шаге производится вывод результатов, то вычисляются необходимые для этого вывода дополнительные переменные для этого момента времени К. После того как окончены все вычисления для момента времени К, текущее время увеличивается на DT единиц времени. Величины, которые были вычислены для момента К, считаются относящимися к моменту J; темпы, вычисленные для интервала KL, рассматриваются как величины на интервале JK; и цикл вычислений снова начинается с расчета уровней. Уравнения в языке DYNAMO имеют вид А.К = V (для уровней, вспомогательных и дополнительных переменных) или B.KL = V (для темпов), где V - означает арифметическое выражение. Всего в DYNAMO имеется около 60 возможных видов правых частей уравнений, что позволяет описывать практически любые зависимости между параметрами модели. Условно уравнения DYNAMO можно разделить на следующие группы. 1. Уравнения для определения новых значений уровней. Правая часть этих уравнений зависит от значения определяемого уровня в момент времени J и значений соответствующих темпов на интервале времени JK. Эти уравнения помечаются в программе буквой L. 2. Уравнения для определения вспомогательных переменных (уравнения помечаются буквой А). 3. Уравнения для определения темпов (помечаются буквой R). 4. Уравнения для присваивания начальных значений параметрам модели (помечаются буквой N) и присваивания значений константам, фигурирующим в описании модели (помечаются буквой С). Кроме этих основных групп уравнений, существуют уравнения для занесения в память машины массивов информации (эти уравнения помечаются буквой Т). Обширный набор специальных функций, которые можно использовать в правых частях уравнений DYNAMO, создает для пользователя большое удобство. Набор включает в себя вычисления элементарных функций анализа (ехр(), sin (x) и др.), функций, позволяющих моделировать случайные воздействия на систему, организовывать запаздывания в контурах обратной связи, активировать и останавливать отдельные процессы в модели. Наконец, стандартные процедуры вывода результатов позволяют получать результаты машинного моделирования в виде удобных таблиц и графиков. При переходе от диаграммы причинных связей к программе на языке DYNAMO обычно используется промежуточная ступень - диаграмма потоков. В отличие от диаграммы причинных связей на диаграмме потоков специальными символами указывается роль каждого параметра в программе на языке DYNAMO. Для этого приняты следующие графические изображения переменных DYNAMO: - Обозначение уровня Х - Обозначение темпа Y - Обозначение вспомогательной переменной В - Обозначение табличных данных Z - Обозначение постоянных факторов W, оценка которых дается экспертизами - Обозначение для неконтролируемых источников материальных потоков Материальные потоки изображаются сплошными линиями, а факт использования значения некоторого параметра при вычислении значения другого параметра отражается на диаграмме потоков ориентированной пунктирной линией, соединяющей изображение этих параметров. Примеры построения диаграммы потоков и программы на языке DYNAMO будут приведены в последующих параграфах. Описание модели на выбранном языке Итак, построена диаграмма причинных связей и выбран язык моделирования. Теперь можно приступить к описанию модели. Этот процесс, однако, не является формальным или по крайней мере чисто формальным. Ведь нам предстоит от диаграммы причинных связей, указывающей лишь на наличие связей между факторами, перейти к установлению точных количественных зависимостей между соответствующими параметрами. Каким же образом устанавливаются эти зависимости? Здесь снова необходимо воспользоваться теми же средствами, которые мы применяли на первом этапе - этапе концептуализации, но трудностей здесь будет больше. Гораздо легче установить экспериментальным путем факт, что изменение одного параметра приводит к изменению другого, чем получить достаточно точную численную характеристику этой зависимости. То же можно сказать и про данные экспертиз. Даже неспециалист может сказать, что число животных должно зависеть от количества и качества корма. Однако для установления точного вида этой зависимости требуется немало труда высококвалифицированных экспертов, обработка большого количества статистических данных и тщательное отсеивание случайных и побочных эффектов. В связи с возникающими в процессе формализации трудностями такого характера особое значение приобретает использование правдоподобных разумных гипотез, которые непременно присутствуют во всякой модели сложной реальной системы. Проверка этих гипотез, выяснение их достоверности сравнением поведения модели с поведением реальной системы зачастую и является основной задачей моделирования. Впрочем, заметим, что такая проверка необходима даже в том случае, если, на наш взгляд, в модель заложены положения, которые нам кажутся бесспорными. Совпадение результатов моделирования с поведением реальной системы не является, конечно, гарантией адекватности модели реальной системе, но тем не менее оно повышает веру пользователя в модель. Хорошие модели реалистичны в деталях и показывают реалистичное поведение. Вернемся к Кэйбабскому плато и рассмотрим диаграмму потоков, построенную в символах, описанных в предыдущем параграфе (рис. 4) На этой диаграмме изображены два уровня: DP - численность оленей на плато и F - количество пищи. Выбор именно этих переменных в качестве уровней определяется уточненной постановкой задачи. Мы будем предполагать, что число хищников в системе - фактор управляемый извне, т.е. будем считать, что количество волков не зависит от поголовья оленей,, а изменяется лишь в результате воздействия внешних факторов (например, охоты на них). Рис. 4. Диаграмма потоков для модели Кэйбабского плато. Отметим, что при построении диаграммы потоков введены дополнительные, по сравнению с диаграммой причинных связей, факторы: FNPD - количество пищи, приходящееся на одного оленя в год, при котором прирост поголовья оленей равен нулю; FR - отношение количества пищи, реально приходящейся на одного оленя, е величине FNPD; GRF - ежегодный процент прироста оленьего стада при данном значении FR. Зависимость GRF от FR имеет нелинейный характер. На рис. 5 приведен график, полученный на основе обработки экспертных оценок и статистического материала. Как видно из этого графика, максимальный прирост поголовья оленей (20 % в год) достигается, если количество пищи, приходящееся на одного оленя, превышает FNPD в полтора раза. Дальнейшее увеличение количества пищи уже ничего не дает для прироста стада. Уменьшение количества пищи ведет к вымиранию оленьего стада, причем максимальная скорость вымирания достигает 30% в год. Рис. 5. Зависимость коэффициента интенсивности естественного прироста от кормового отношения. В модели имеются еще три переменные величины, изменения которых также имеют нелинейных характер и задаются соответствующими таблицами. На рис. 6 приведена зависимость количества оленей, уничтожаемых одним хищником в год (величина DKR), от средней плотности популяции оленей на плато (величины DD). При плотности оленей 0,025 на акр, т.е. один олень на 16 га, величина DKR достигает своего максимума (60 оленей в год). Рис. 6. Зависимость числа оленей, уничтожаемых одним хищником за год, от плотности оленей. На рис. 7 дается график зависимости количества пищи, потребляемой одним оленем в год FCPD, от величины FR. Из рисунка видно, что в случае нехватки пищи она полностью съедается оленями. Когда же образуется избыток пищи, по сравнению с величиной FNPD, поедание кормов будет неполным, что благотворно сказывается на восстановлении пастбищ на плато. Рис. 7. Зависимость корма, потребляемого одним оленем, от кормового отношения. Наконец, последний график, представленный на рис. 8, показывает изменение времени восстановления пастбищ FRT от отношения реального количества пищи на пастбищах F к максимально возможному количеству пищи FN. Рис. 8. Зависимость времени воспроизводства корма от отношения имеющегося корма к максимально возможному количеству корма. Смысл параметра "время восстановления пастбища" (FRT) нуждается в дополнительном пояснении. Пусть в некоторый момент времени количество пищи на пастбище меньше максимально возможного количества FN на величину дельта F. Тогда FRT представляет собой то время, за которое эта разница уменьшается в 2,73 раза (точнее, в е раз, где е - основание натуральных логарифмов). Выбор такой величины в качестве одного из параметров модели объясняется тем, что экспериментальные кривые восстановления травяного покрова за достаточно короткий промежуток времени, не превышающий одного сезона, аппроксимируются экспоненциальной зависимостью от времени. Величина FRT аналогична постоянной времени, встречающейся при решении линейных дифференциальных уравнений первого порядка. Следующий шаг реализации модели заключается в выборе языка и написании программы для ЭВМ. Составление программы - это завершение точной формализованной постановки задачи и выбора метода ее решения. Все зависимости, от которых до этого еще можно было говорить, не конкретизируя их вида, теперь должны быть описаны в виде формул на соответствующем языке программирования. При программировании задаются значения всех постоянных величин, включаемых в модель, точно задаются данные, описывающие начальное состояние модели, определяется, что и в каком виде нужно получить с выводных устройств ЭВМ как результат моделирования. Ниже приведена программа модели Кэйбабского плато на языке DYNAMO. Программа состоит из четырех частей (названия этих частей помещены в качестве замечаний рядом со словами NOTE). NOTE <*> K A I B A B P L A T E A U M O D E L <*> NOTE NOTE ***** D E E R P O P U L A T I O N ***** NOTE L DP.K = DP.J + DT*(DNGR.JK-DPR.JK) N DP = DPI C DPI = 4000 R DNGR.KL = DP.K * GRF.K A GRF.K = TABHL (GRFT, FR.K, 0, 1.5, .25) A FR.K = FPD.K / FNPD C FNPD = 1 T GRFT = -.3/-.27/-.22/-.15/0/.15/.2 A FPD.K = F.K / DP.K R DPR.KL = DKR.K * PP.K A DKR.K = TABHL (DKRT, DD.K, 0, .025, .005) T DKRT = 0/3/13/32/51/56 A DD.K = DP.K/AREA C AREA = 300000 NOTE ***** NOTE ***** P R E D A T O R P O P U L A T I O N ***** NOTE ***** A aa.K = PPI*CLIP(EXP(-GC*(TIME.K-BTIME)),1,TIME.K,BTIME) A PP.K = CLIP(100,aa.K,TIME.K,1950) C PPI = 266 N PP = PPI C GC = .2 C BTIME = 1904 NOTE ***** NOTE ***** F O O D S E C T O R ******* NOTE ***** L F.K = F.J + DT * (GR.JK - FC.JK) N F = FI C FI = 350000 R GR.KL = (FN - F.K) / FRT.K C FN = 350000 A FRT.K = TABLE (FRTT, F.K / FN, 0 , 1, .25) T FRTT = 14/6/3/1.5/1 R FC.KL = DP.K * FCPD.K A FCPD.K = TABHL (FCPDT, FR.K, 0, 1.5, .25) T FCPDT = 0/.25/.5/.75/1/1.12/1.2 NOTE ***** NOTE ***** C O N T R O L C A R D S **** NOTE ***** N TIME = TIMEI C TIMEI = 1880 SPEC DT = .1/LENGTH=1980/PRTPER=5 PRINT TIME/DP/F/PP PLOT DP/F/PP INPUT TIMEI/LENGTH/PPI/DPI INTAB DKRT(6) Первая часть программы содержит уравнения, описывающие изменение численности оленей, вычисление соответствующих вспомогательных переменных, занесение в память необходимых для работы этой части программы таблиц, присвоение переменным начальных значений, присвоение значений константам. Поясним уравнения этой части L DP.K = DP.J + (DT)(DNGR.JK - DPR.JK) Символ L ( начальная буква английского level - уровень) указывает, что данное уравнение - уравнение уровня. Количество оленей в текущий момент времени - уровень DP.K - равняется значению этого уровня в предыдущий момент времени (DP.J) плюс разность между темпом увеличения поголовья оленей за счет естественного прироста на интервале JK (величина DNGR.JK) и темпом убыли поголовья за счет уничтожения хищниками на том же интервале JK ( величина DNGR.JK) и темпом убыли поголовья за счет уничтожения хищниками на том же интервале JK (величина DPR.JK), умноженная на величину этого временного интервала DT. Следующее уравнение N DP = DP1 осуществляет присвоение уровню DP начального значения, равного константе DP1. Символ N указывает на операцию присвоения начального значения. Третье уравнение присваивает константе DP1 значение 4000. Четвертое уравнение программы R DNGR.KL = DP.K * GRF.K указывает, что численный прирост поголовья на следующем (за данным) интервале времени KL равен числу оленей в настоящий момент DP.K, умноженному на коэффициент естественного прироста. Символ R отмечает, что это уравнение - уравнение темпа. Восьмое уравнение программы T GRFT = -.3/-.27/-.22/-.15/0/.15/.2 показывает, что в памяти машины организовывается одновходовая таблица некоторой функции (табличные значения перечислены в правой части уравнения). Идентификатор GRFT является названием таблицы. В пятом уравнении программы A GRF.K = TABHL (GRFT, FR.K, 0, 1.5, .25) эта таблица используется для присвоения значения величине GRF.K - коэффициенту естественного прироста оленей. Смысл аргументов функции TABHL в правой части уравнения следующий: GRFT - название соответствующей таблицы, FR.K - значение входа в таблицу, т.е. то значение, от которого нам нужно определить величину функции GRF. Следующие два аргумента (0 и 1,5) - соответственно нижняя и верхняя границы области изменения независимой переменной. Последний аргумент (0,25) указывает величину шага таблицы. Если значение входа не совпадает ни с одним из узлов таблицы, то значение функции получается с помощью линейной интерполяции. Аргументам, выходящим за пределы таблицы, присваиваются соответствующие крайние значения. Вторая часть описывает численность популяции хищников. Поясним первое уравнение этого раздела программы. A PP.K = PP1 * CLIP (EXP(-GC * TIME.K - BTIME)),1, TIME.K, BTIME) Выражение PP.K, стоящее в левой части уравнения, обозначает количество хищников в текущий момент времени. В правой части уравнения переменная TIME.K обозначает номер текущего момента времени, в данном случае года. Функция EXP (A * (B - C)) в обычной математической записи означает е в соответствующей степени. Специальная функция CLIP определяется следующей формулой: { а, если с >= d, CLIP (a,b,c,d) = { b, если c < d. Таким образом, предполагается, что численность до момента времени BTIME - постоянна и равна PP1, а затем экспоненциально убывает со скоростью, определяемой показателем GC. Третья часть программы содержит уравнения, отражающие изменения количества пищи на плато. Наконец, четвертый раздел программы - служебный. В нем (стандартным для DYNAMO образом) указано начальное значение времени (TIME1), величина временного шага (DT), конец интервала времени, в течение которого моделируется явление (LENGTH), шаг выдачи результатов на графики (PLTPER). Последняя строка этого раздела указывает, какие величины выводятся на графики, какими символами отмечаются точки этих кривых. Масштабы по оси ординат для всех графиков, кроме графика величины PP, выбираются автоматически. В результате вычислений, проведенных авторами этой модели на ЭВМ, были получены результаты, которые мы сейчас обсудим. Все моделирование проводилось с двумя целями. Первая цель - показать, что построенная модель хорошо отражает реальное изменение численности оленьего стада на Кэйбабском плато. Вторая цель - проверить действие различных факторов, влияющих на изменение поголовья оленей. Рассмотрим несколько упрощенных вариантов модели. В первом варианте предполагается, что количество корма не ограничено. Это приводит к тому, что естественный прирост поголовья оленей достигает своей максимальной величины - 20% в год. Изменение численности оленей зависит в этом варианте только от количества хищников. При этом предположении диаграмма причинных зависимостей сводится к диаграмме, приведенной на рис. 9. Рис. 9. Модель Кэйбабского плато. Диаграмма причинных связей в предположении о неограниченном количестве корма. На рис. 10 приведена зависимость количества оленей при постоянном числе хищников (50 особей). На этом же рисунке приведены динамика изменения количества оленей, уничтожаемых одним хищником в год. Рис. 10. Изменение численности оленей со временем (50 волков, корм не ограничен. Рис. 11 показывает зависимость численности популяции оленей от времени, когда на плато присутствуют 300 хищников. Графики показывают, что в первом случае популяция оленей неограниченно растет, а во втором - убывает. Как показали эксперименты в предположении, что количество пищи неограниченно, равновесие системы достигается в случае, если число хищников равно 26. При этом число оленей составляет около 4600 голов. Следующий упрощенный вариант модели подразумевает постоянство количества пищи. Качественное поведение системы незначительно отличается от предыдущего случая. Рис. 11. Изменение численности оленей со временем (300 волков, корм неограничен). На рис. 12 приведены результаты имитации процесса восстановления кормовой базы на плато в предположении, что число оленей постоянно (4000 голов). Рис. 12. Эффект восстановления кормовой базы при сохранении постоянной численности оленей (4000 голов) в 1904 г. Рис. 13 показывает, как влияет на систему отстрел волков, начавшийся примерно в 1904 г. При этом предполагается, что время воспроизводства корма постоянно и равно пяти годам. Мы видим, что количество хищников сокращается, как это предполагалось соответствующим уравнением программы. Оно приводит к резкому возрастанию численности оленей, в результате чего количество пищи столь же резко падает. Количество оленей при этом несколько уменьшается, после чего снова достигается равновесие. Рис. 13. Эффект массовой охоты на волков (время воспроизводства корма постоянно - 5 лет). Рис. 14. Эффект массовой охоты на волков в 1904 г. (без предположения о постоянстве времени воспроизводства корма). Следующий рисунок (рис. 14) представляет результаты моделирования несколько более реальной ситуации - снимается предположение о постоянстве времени восстановления корма. Качественно кривые выглядят так же, как и в предыдущем случае, однако равновесные значения численности поголовья оленей и количества пищи уменьшились вдвое. Интересно сравнить этот рисунок с рис. 2, показывающим изменение реальной численности оленей на Кэйбабском плато в период с 1905 по 1939 г. после почти полного отстрела хищников. Этот график приведен американским ученым Корманди в книге "Концепции экологии". На нем отмечены те моменты (1910, 1920, 1923 гг.), когда специалисты выражали беспокойство по поводу будущей судьбы чрезвычайно быстро растущего оленьего стада на плато. Действительно, ход кривой показывает, что в 1925 г. наступил кризис - за две зимы численность оленей сократилась на 60%. Сильно подорванная кормовая база плато так и не смогла восстановиться до нормального уровня, при котором на нем могло, по оценкам специалистов, прокормиться около 30 тыс. оленей. В результате численность оленей продолжала убывать; с 40 тыс. в 1925 г. упала к 1939 г. до 10 тыс. А между тем, как показали расчеты на модели, если бы в 1922 г. с помощью отстрела оленей численность стада была бы доведена до 25 тыс., то это поголовье могло бы находиться в равновесии долгие годы (рис. 15). Если бы было оставлено даже 50 тыс. оленей, то все равно часть бы вымерла от недостатка пищи, однако тогда кормовая база смогла бы восстановиться и вновь было бы достигнуто равновесие при двадцатипятитысячном стаде (рис. 16). Кто знает, если бы в 1922 г. были бы проведены эти расчеты и приняты соответствующие меры, то может быть Кэйбабское стадо оленей не пришло бы к плачевному состоянию к 1939 г. Интересно отметить, что экологи предсказывали нечто подобное и близкое по цифрам (см. рис. 2). Рис. 15. Влияние регулирования численности оленей на уровне 25 тысяч. Рис. 16. Влияние регулирования численности стада на уровне 50 тысяч. Рис. 17 демонстрирует эффект уменьшения времени восстановления кормовой базы: численность оленей к 1950 г. возрастает вдвое, количество пищи также увеличивается и ее хватает для стада возросшей численности. Рис. 17. Эффект уменьшения времени восстановления корма в 1950 г. Влияние введения в систему сотни хищников изображено на рис. 18. Результаты моделирования показывают, что хотя при этом количество оленей уменьшается (опять примерно вдвое), количество пищи резко увеличивается, достигая за 10 лет своего максимально возможного уровня. Рис. 18. Результат возвращения в систему 100 волков (1950 г.). Конечно, описанная модель не учитывает много факторов, могущих повлиять на динамику системы. Так, например, не учитывается то обстоятельство, что хищники, как правило, уничтожают более слабых животных, участие которых в репродуктивной деятельности стада нежелательно с биологической точки зрения. Следует отметить, что высокий уровень агрегирования модели, при котором предполагается, что на плато, поросшем "осредненной" травой, пасется стадо из одинаковых "средних" оленей, приводит к ошибкам такого рода: в модели "средние" олени в тяжелых условиях вымирают поголовно, в жизни же почти всегда находятся особи, которые в силу либо свои личных качеств, либо в силу каких-то локальных особенностей области обитания переживают трудный период и обеспечивают затем сохранение вида. Более того, в генетике предполагается, что стабильный процент мутантов в популяциях земных видов как раз и обеспечивает "поставку" особей, которые при резком неблагоприятном изменении условий жизни могут оказаться жизнеспособней основной массы потомства. Описанные результаты тем не менее показывают, что из модели можно получить много ценной информации как для науки, так и для практической деятельности по охране природы. Отдельные закономерности и приемы, использованные в этой модели, относятся не только к оленям и волкам, а могут быть с успехом применены для создания моделей, с помощью которых можно решать ряд сложных экологических проблем. Часть II. Описание языка моделирования DYNAMO Глава I. Основные понятия 1.1 Понятие уровня Основным понятием, используемым в языке DYNAMO, является понятие уровня. Например, рассмотрим количество жидкости в баке как функцию от скорости поступления жидкости и времени. Если в баке час назад было 100 л и поток составляет 1л/мин., то сейчас в баке будет 160 л жидкости. В виде уравнения это можно записать так: ┌ количество ┐ ┌ количество ┐ ┌ прошедшее┐ ┌ скорость ┐ │ жидкости в │ = │ жидкости в │ + │ время │ * │ поступления │ └ баке сейчас┘ └ баке раньше┘ └ ┘ └жидкости в бак┘ Настоящий момент времени в DYNAMO принято обозначать индексом К, прошлый - индексом J, будущий - L. Время, прошедшее между J и K - за DT. Если обозначить количество жидкости в баке за Х, то можно переписать уравнение так: ┌скорость поступления ┐ X.K = X.J + DT * │жидкости в бак │ └ ┘ (Десятичная точка отделяет временной индекс от идентификатора X). Этот тип уравнений называется уровневым, а величина в его левой части (X) - уровнем. В этом примере скорость потока жидкости в бак считается постоянной, в то время как на практике она переменна. Поэтому, чтобы получить правильное решение, надо временной интервал разбить на некоторое число меньших интервалов, в пределах каждого из которых значение скорости можно считать постоянным. После этого повторить вычисления для каждого интервала (DT следует положить равным длине интервала). Существует простой метод аппроксимации функции скорости, дающий хорошие результаты если не требуется большая точность вычислений (метод Эйлера). С формальной точки зрения язык DYNAMO как раз и описывает процесс численного интегрирования системы дифференциальных уравнений методом Эйлера. Точность вычислений регулируется соответствующим подбором величины DT. 1.2 Консервативные системы. Понятие скорости (темпа) Пример с потоком жидкости в бак является типичным для многих ситуаций, где что-то перемещается с места на место без потерь и увеличения во время процесса. Подобными примерами являются поток товаров на склад, люди в потоке рабочей силы, ток в конденсаторе. Такой тип систем назовем консервативными подсистемами, а потоки - темпами. В консервативных системах темпы изменения значений уровней имеют вид суммы или разности нескольких темпов. В тех случаях, когда темп только увеличивает значение уровня и никогда его не уменьшает, он трактуется как приток из внешнего источника. Если же темп только уменьшает значение уровня, то его рассматривают как поток к внешнему стоку. Темпы вычисляются в текущий момент времени К для интервала KL из алгебраических уравнений, включающих значения переменных в тот же самый момент времени (K). Пример уравнения темпа: R.KL = LEV.K/DEL 1.3 Неконсервативные системы. Понятие дополнения Неконсервативные называются все системы, не являющиеся консервативными. Здесь интегралы также вычисляются в уровневых уравнениях, но темпы изменения значений уравнений носят более сложный характер, чем сумма нескольких темпов (как это имеет место в консервативных системах). Типичным примером может быть следующее уравнение: LEV.K = LEV.J + DT * (R.JK - LEV.J)/CONST Здесь скорость изменения значения уровня LEV есть разность между величинами с разными индексами, деленная на константу. Такая запись уравнений вполне допустима. Алгебраические связи в неконсервативной системе вычисляются с помощью уравнений дополнения. Значения этих уравнений вычисляются в момент времени К с использованием значений уровней и других дополнений в момент времени К. Уравнения дополнения вычисляются так, чтобы ни одно из них не использовалось в другом уравнении раньше, чем будет получено его значение. Дополнения часто используются как составные части уравнений темпов. 1.4 Последовательность вычислений Первыми в текущий момент времени К вычисляются уровни. Они зависят от собственного значения в предыдущий момент времени J, а так же от дополнений, вычисленных в момент времени J, и скоростей, вычисленных на интервале JK. Так как все величины для момента времени J и интервала JK вычислены, то можно вычислить значения уровней в момент времени К. Затем вычисляются дополнения для момента времени К с использованием значения уже вычисленных в момент времени К уровней и дополнений. Последними для интервала KL вычисляются значения скоростей с использованием значений уровней и дополнений, вычисленных в момент времени К. Порядок вычисления скоростей несущественен, так как во всех уравнениях используются их значения в интервале JK. Таковы общие правила. Данный компилятор распознает любые зависимости между величинами, входящими в уравнения (в том числе и более сложные, чем описаны выше), если только между ними нет перекрестной зависимости, когда одна величина определяется через другую а та в свою очередь через нее и т.п. Когда значения скоростей вычислены время увеличивается на DT; все величины, вычисленные в момент времени К и на интервале KL, рассматриваются как вычисленные в момент J и на интервале JK соответственно. Затем начинается новая итерация вычислений. 1.5 Запись уравнений Основной формат записи уравнения: <тип><пробелы><величина> = <выражение> <тип> - это буква, определяющая тип уравнения: L - уравнение уровня, A - уравнение дополнения, R - уравнение скорости (темпа) N - уравнение задания начальных условий, C - уравнения задания констант, T - уравнения задания таблиц. Буква, указывающая тип уравнения должна быть первой в строке! <величина> - идентификатор с временным индексом, отделенным от него десятичной точкой. В левой части уравнений могут быть только индексы K и KL(в уравнениях темпов).В уравнениях N C и T типов временные индексы не нужны.Идентификатор может состоять из букв и цифр, причем первый символ - буква, идентификатор может иметь любую длину, но только первые десять символов распознаются компилятором. <выражение>- комбинация множителей и термов, включающих пере- менные, числа и обращения к функциям. В DYNAMO ис- пользуются операции сложения, вычитания, умноже- ния, деления, обозначаемые символами + - * /. Порядок выполнения действий слева направо: сначала выполняются действия в скобках, производятся обра- щения к подпрограммам- функциям, выполняются в по- рядке следования операции умножения и деления, а затем сложения и вычитания. В DYNAMO имеются следующие встроенные функции: 1. Тригонометрические : COS, SIN, EXP, LOGN, SQRT. 2. Выбор значения : CLIP,MAX,MIN, SWITCH,TABLE,TABLH. 3. Временные : PULSE,RAMP,SAMPLE,STEP. 4. Задержки : SMOOTH, DELAY3, DLINF1, DLINF3. 5. Датчики случайных чисел : NOISE,NORMRN. Аргументами функций могут быть любые выражения,а также другие функции. Имена функций должны писаться большими латинскими буквами. 1.6 Задание начальных условий Чтобы запустить цикл вычислений, надо задать начальные значения уровней. Они задаются при помощи уравнений задания начальных условий. Эти уравнения имеют тип N и составляются по общим для всех уравнений правилам, только величины входят в них без временных знаков. Например, N ABC = 1E-2 N KLM = RRR * RK Во втором случае значение величины KLM задается при помощи величин RRR и RK. Если это дополнения или скорости, то для них тоже могут быть заданы начальные значения. Если же значение их не задано, то будут использоваться описывающие их уравнения как N-уравнения. Этот процесс продолжается до тех пор пока не будут найдены начальные значения величин как конечные уравнения. Для задания констант - параметров модели используются C-уравнения. Пример: C WEIGHT = 1.2E10 ГЛАВА 2 Описание языка 2.1 Алфавит языка 1. Заглавные и прописные буквы латинского алфавита 'A B C ... Z a b c ... z' Заглавные и прописные буквы в идентификаторах различаются (рекомендуется использовать в программах только заглавные латинские буквы). 2. Цифры '0 1 2 3 4 5 6 7 8 9' 3. Знаки арифметических операций и скобки '+ - * / = ( )' 4. Десятичная точка '.' 5. Пробел. 2.2 Грамматика языка Ниже приводится грамматика языка в форме Бэкуса - Наура. Нетерминалы обозначены < > . SPEC L A R T N C B PRINT PLOT INPUT INTAB EXTRN /* заpезеpвиpованные слова */ <Программа> : <строки> <строки> : : <строка><строки> ; <строка> : '\n' /* переход на новую строку */ | <уравнение> '\n' | <спецификация> '\n' | <вывод> '\n' | <ввод параметров> '\n' | <внешние переменные > '\n' ; <уравнение> : | | | | | ; : L <переменная> '=' <выражение>; : R <переменная> '=' <выражение>; : A <переменная> '=' <выражение>; : N <переменная> '=' <выражение>; : C <идентификатор константы> '=' <константа>; : T <идентификатор таблицы> '=' <список чисел>; <список чисел> : <число> | <число> '/' Nl <список чисел> ; <выражение> : <выражение> '+' <выражение> | <выражение> '-' <выражение> | <выражение> '*' <выражение> | <выражение> '/' <выражение> | '(' <выражение> ')' | '-' <выражение> /* унарный минус */ | <терм> ; <терм> : <имя функции> '(' <список аргументов> ')' | <имя функции> | <переменная> | <число> | <целое> ; <список аргументов> : <выражение> | <выражение> ',' <список аргументов> ; <константа> : <знак> <число> | <знак> <целое> ; <переменная> : <идентификатор> | <идентификатор> <время> ; <время> : '.J' | '.K' | '.KL' | '.JK' ; <знак> : '+' | '-' | '*' | '/' ; <спецификация> : SPEC Nl <список спец> <список спец> : <список1> | <список1> '/' Nl <список спец> ; <список1> : <переменная> '=' <константа> ; <вывод> : PRINT <список на печать> | PLOT <список на график> ; <список на график> : <переменная> | <переменная> '/' Nl <список на график> ; <список на печать> : <пер_на_печать> | <пер_на_печать> '/' Nl <список на печать> ; <пер_на_печать> : <переменная> <поле> ; <поле> : '(' <целое> ',' <целое> ')' | /* пусто */ ; <ввод параметров> : INPUT <список параметров> | INTAB <список таблиц> ; <список параметров> : <переменная> | <переменная> '/' Nl <список параметров> ; <список таблиц> : <переменная> '('[целое]')' | <переменная> '('[целое]')''/'Nl <список таблиц> ; Nl : | '\n' ; <число> : <целое> |<целое>'.'<цел> |<целое>'.'<цел>'E'<цел> |<знак>'.'<цел>; <целое> : <знак><цифра> | <знак><целое><цифра> ; <цел> :<цифра>|<цел><цифра> <цифра> :0|1|2|3|4|5|6|7|8|9; <идентификатор> : <буква>| <идентификатор><буква> | | <идентификатор><цифра>; <буква> :A|B|C|D|E|F|G|H|I|J|K|L|M|N|O|P|R|S|T|U|V|W|X|Y|Z| a|b|c|d|e|f|g|h|i|j|k|l|m|n|o|p|r|s|t|u|v|w|x|y|z; Комментарий начинается с новой строки ( с первой колонки) служебным словом NOTE. Эта строка компилятором игнорируется. Пример программы на языке DYNAMO, написанный в соответствии с настоящей грамматикой, приводится в первой части книги. 2.3 Идентификатор переменной Идентификатор переменной может состоять из любого количества символов, но компилятор учитывает только первые 10. Идентификаторы DT, LENGTH, TIME, PRTPER и имена встроенных функций зарезервированы в DYNAMO и не могут применяться пользователем для обозначения его собственных величин. 2.4 Числа. Число может состоять не более чем из 15 цифр (мантисса ). Порядок от E-308 до E+308. Примеры: 18 -223 1.2464E-15 .234 2.5 Временные индексы Уровни и дополнения имеют временные индексы из одной буквы J или K, скорости же имеют индекс из двух букв JK или KL. Текущему моменту времени соответствует индекс K прошедшему J. Вычисляемой величине (в левой части уравнения) всегда соответствует индекс K или KL в зависимости от типа уравнения. В уравнениях задания начальных условий (тип N) все величины пишутся без индексов. Уравнение C и T типов подробно описаны в грамматике. В таблице приведены индексы, используемые в различных уравне- ниях DYNAMO. ┌───────────────┬───────────┬────────────────────────────────────┐ │ Тип уравнения │ Временный │ Временный индекс величины в правой │ │ │ │ части уравнения │ │ │ индекс ├──────┬──────┬─────┬─────┬─────┬────┤ │ │ слева │ L │ A │ R │ C │ N │ T │ ├───────────────┼───────────┼──────┼──────┼─────┼─────┼─────┼────┤ │ L уровень │ K │ J │ J │ JK │ - │ - │ - │ │ А дополнение │ K │ J K │ J K │ JK │ - │ - │ - │ │ R темп │ KL │ J K │ J K │ JK │ - │ - │ - │ │ C константа │ - │ - │ - │ - │ - │ - │ - │ │ N нач. условие│ - │ - │ - │ - │ - │ - │ - │ │ T ввод таблицы│ - │ - │ - │ - │ - │ - │ - │ └───────────────┴───────────┴──────┴──────┴─────┴─────┴─────┴────┘ Следует помнить что каждая величина в правой части уравнения должна быть описана в левой части какого-либо уравнения. Имеются два исключения из этого правила: величины DT и TIME. Величина DT определяется в спецификации и может быть использована как константа в любом выражении. Величина TIME встроена в язык и нигде не определяется в программе (кроме случая когда необходимо начать счет не с нулевого момента времени, тогда начальное значение TIME присваивается с помощью N уравнения). Она может быть использована в качестве уровня в правой части любого уравнения. 2.6 Задание начальных условий Все L-переменные должны иметь начальные условия, задаваемые c помощью N-уравнений. Дополнения и скорости, употребляемые в N-уравнениях, могут иметь начальные значения, определяемые пользователем или по умолчанию. 2.7 Спецификации запуска Существует 3 величины, определяющие важные параметры запуска модели. DT - интервал времени TIME между моментами J и K (шаг интегрирования) LENGTH - значение времени TIME, когда следует остановить запуск. PRTPER - интервал времени, через который следует печатать результат. Существует верхний предел на величину DT, определяющийся величиной задержек первого порядка, действующих в системе. DT должна быть вдвое меньше минимальной величины задержки первого порядка, существующих в модели. Для вывода на печать необходима спецификация вида: PRINT TIME/DP(12,4)/PP Эта спецификация указывает компилятору, что следует вывести на печать ( в файл с расширением .RES) величины TIME, DP и PP, причем у величин TIME и PP будут выведены по 8 знаков и 2 знака после запятой, (по умолчанию),а у величины DP - 12 знаков и 4 знака после запятой - на это указывает спецификация DP(12,4). Для графического вывода необходима спецификация вида: PLOT DP/PP. Эта спецификация указывает компилятору, что величины DP и PP следует вывести на график. В спецификации PLOT может быть не более 5 величин Если будет указано больше, то выведены будут только первые 5. 2.8 Встроенные функции Секущая функция CLIP (P,Q,R,S) ┌ P если R >= S CLIP = │ └ Q если R < S Задержки SMOOTH - задержка первого порядка в материальном потоке. Пример употребления A SMTHX.K = SMOOTH(K.JK, SMTM), что является сокращенной записью следующих уравнений: L $L1.K = $L1.J + DT * (X.JK - $R1.JK) N $L1 = X * SMTM A SMTHX.K = $L1.K / SMTM R $R1.KL = $L1.K / SMTM Уровень $L1 и темп $R1 продуцируются транслятором автоматически, запаздывание есть константа. DELAY3 - задержка третьего порядка в материальном потоке. Пример употребления R R.KL = DELAY3(Q.JK, DEL) что является сокращенной записью следующей цепочки уравнений: L $L1.K = $L1.J + DT * (Q.JK - $R1.JK) N $L1.K = Q * ( DEL/3 ) R $R1.KL = $L1.K / (DEL/3) L $L2.K = $L2.J + DT * ($R1.JK - $R2.JK) N $L2.K = Q * ( DEL/3 ) R $R2.KL = $L2.K / (DEL/3) L $L3.K = $L3.J + DT * ($R2.JK - R.JK) N $L3.K = Q * ( DEL/3 ) R R.KL = $L3.K / (DEL/3) DLINF1 - задержка первого порядка в информационном потоке Пример употребления A RVX.K = DLINF1( X.K, TRX) что является сокращенной записью уравнений: R $R1.KL = (X.K - $L1.K) / TRX L $L1.K = $L1.J + DT * $R1.JK N $L1 = X A RVX.K = $L1.K DLINF3 - задержка третьего порядка в информационном потоке Пример употребления A RVX.K = DLINF3(X.K, TRX) что является сокращенной записью уравнений: R $R1.KL = (X.K - $L1.K) / (TRX/3) L $L1.K = $L1.J + DT * $R1.JK N $L1 = X R $R2.KL = ($L1.K - $L2.K) / (TRX/3) L $L2.K = $L2.J + DT * $R2.JK N $L2 = X R $R3.KL = ($L2.K - $L3.K) / (TRX/3) L $L3.K = $L3.J + DT * $R3.JK N $L3 = X A RVX = $L3.K MAX(P,Q) Функция максимума ┌ P если P>=Q MAX (P,Q) = │ └ Q если P=Q Случайные числа NOISE -дает случайные числа равномерно распределенные в интервале: [-1/2 ,+1/2] NORMRN(Q,R) - нормально распределенные числа с математическим ожиданием Q и среднеквадратическим отклонением R Источник импульсов. PULSE (P,Q,R) - серия импульсов шириной DT и высотой P, первый импульс появляется в момент Q, последующие с интервалом R. Функция наклона RAMP ( SLP,STRT) ┌ 0 если TIME <= STRT RAMP =│ └ сумме всех SLP*DT от STRT до TIME.K Функция выборки SAMPLE(P,Q,R) равняется P в моменты выборок с интервалом Q и сохраняет свое значение до следующего момента выборки через интервал Q, R - начальное значение SAMPLE до первого выбора. Шаговая функция ┌ 0 если TIME < Q STEP (P,Q) =│ └ P если TIME >=Q Функция переключатель ┌ P , если R=0 SWITCH (P,Q,R) = │ └ Q , если R<>0 Функции ввода таблиц V = TABLE(NAME,P,N1,N2,N3) - с помощью этого уравнения переменной V присваивается значение F(P). Функция F(P) задается таблицей NAME; N1, N2 - соответственно нижняя и верхняя границы области изменения независимой переменной. N3 - шаг таблицы. P - независимая переменная. Если значение P не совпадает с одним из узлов таблицы, то F(P) получается с помощью линейной интерполяции. В этом уравнении аргумент P обязательно должен удовлетворять неравенству N1<=P<=N2 V = TABHL(NAME,P,N1,N2,N3) - это уравнение отличается от предыдущего тем, что аргумент P может принимать произвольное значение.Если PN2, то V принимает значение F(N2). Г Л А В А 3. Инструкция пользователя. Требования данного компилятора к конфигурации ЭВМ: - IBM PC, PC /XT, PC/AT или совместимые; - DOS 2.0 и выше; - Любой графический адаптер; Не менее 512 Kб RAM. 3.1 Комплект поставки (версия DYNAMO BASED EDITION) DYN.EXE - программа, обеспечивающая трансляцию DYNAMO-программ в формат .EXE DYNINST.EXE - установщик для DYN.EXE DYNAMO.H - файл с объявлениями для C-компилятора DYNAMO.DOC - этот файл SAMPLES.ARJ - примеры DYNAMO-программ INCLUDE.ARJ - файлы заголовков для C-компилятора QCL.ARJ - библиотеки QCSUP.ARJ - прочие Quick C файлы ARJ.EXE - архиватор; LIB.ARJ - библиотеки для DYNAMO 3.2 Запуск транслятора Использование транслятора DYNAMO-программ: DYN.EXE 3.3 Выполнение программы, построенной с помощью DYNAMO транслятора В данной версии транслятора возможно тестирование модели без перекомпилирования программы. Для этого необходимо где либо в программе указать спецификацию INPUT или INTAB. Например: INPUT DP/PP INTAB DKRT(6) Эти спецификации укажут транслятору, что необходимо предусмотреть возможность динамического изменения величин DP и PP и таблицы DKRT (6 - размер таблицы). Если в программе отсутствуют спецификации PLOT, INTAB или INPUT, весь результат работы программы - таблица в файле с таким же именем, но с расширением .RES. Примечание! Должна быть хотя бы одна опция PRINT или PLOT в Вашей программе, иначе Ваша программа не будет производить никакого вывода. Полученную таблицу можно посмотреть в любом редакторе. Если в вашей программе есть спецификация PLOT, то после выполнения счета на экране появится меню вида: Common Chart Tour of Charts (1) Quit Выбрав в нем опцию Common Chart Вы сможете посмотреть все величины в спецификации PLOT на одном графике в относительных измерениях. Выбрав Tour of Charts вы сможете просмотреть все величины последовательно на графиках в абсолютных измерениях. Quit - выход. Если в программе присутствуют спецификации ввода (INPUT или INTAB), то сначала появится меню Execution New Parameters (2) New Tables Exit Execution - выполнение с данными значениями параметров. New Parameters и New Tables - ввод-изменение значений параметров и таблиц. После изменения значений снова появляется это же меню. Exit -конец счета модели. В это же меню модель будет возвращаться после каждого сеанса интегрирования с восстановленными значениями всех параметров и таблиц. ВНИМАНИЕ! 1. Не следует указывать среди параметров PLOT величины, являющиеся константой по времени. 2. На графиках не предусмотрено отображение отрицательных величин, они автоматически пересчитываются (сдвигаются в положительную область). Эти недоработки будут устранены в последующих версиях. 3.4 Сообщения об ошибках "Equation table owerflow" - переполнена таблица уравнений. У вас очень большая модель. Чтобы транслятор мог обработать столько уравнений, необходимо изменить значение переменной LEN_EQU_ TABLE в файле на большее. Сейчас оно равно 400 уравнений. "syntax error" - синтаксическая ошибка. "This name already defined as not a level" - эта переменная уже была определена не как уровень. Проверьте "This name already defined as not a rate" - эта переменная уже была определена не как темп. Проверьте. "This name already defined as not an append variable" - эта переменная уже была определена не как вспомогательная переменная. Проверьте. "This name already defined as constant or table" - это имя уже определено как константа или таблица. "This name already defined as not a constant" - это имя уже определено не как константа. "This name already defined as not a table" - это имя определено уже не как таблица. "[first] argument of this function must be variable" - [первый] аргумент этой функции должен быть переменной. "can't open input file" - не могу открыть входной файл. Возможно вы ошиблись в написании имени файла. "can't open output file" - не могу открыть выходной файл. "can't open error file" - не могу открыть файл для сообщений об ошибках. "Syntax errors were found." - были обнаружены синтаксические ошибки. "Invalid usage of DELAY" - неправильное использование функции задержки: следует использовать ее так: R V.KL = DELAY(IN. JK,DEL) "Invalid usage of DELAY3" - неправильное использование функции задержки третьего порядка, следует использовать ее так: R V.KL = DELAY3 (IN.JK,DEL) "symbol table overflow" - переполнена таблица символов. "out of memory" - нет памяти для дерева. "Error:cross-dependence between equations: check equations - фатальная ошибка: перекрестная зависимость между уравнениями. 3.5 Предупреждения. "On the left side of equation must be ".K" time" - в левой части уравнения должна быть переменная с временным суффиксом ".K" "On the left side of equation must be level" - в левой части уравнения должен быть уровень. "This level already defined.New definition ignored" - этот уровень уже определен, новое определение игнорируется. "On the left side of equation must be ".KL" time" - в левой части должен быть временной суффикс ".KL" "On the left side of equation must be rate" - в левой части уравнения должен быть темп. "This rate already defined.New definition ignored" - этот темп уже определен, новое определение игнорируется. "On the left side of equation must be no time " - в левой части не должно быть временного суффикса ( N - уравнение). "This variable already initialised." "New initialisation ignored" - эта переменная уже инициализирована N-уравнением, новая инициализация игнорируется. "On the left side of equation must be name" - в левой части уравнения должно быть имя (C-уравнение) "This constant already defined. New definition ignored" - эта константа уже определена, новое определение игнорируется. "This table already defined. New definition ignored" - эта таблица уже определена, новое определение игнорируется. "incorrect number of arguments" - неверное количество аргументов у функции. "print option alredy specificated" - спецификация PRINT уже сделана, эта игнорируется. "plot option alredy specificated" - спецификация PLOT уже сделана, эта игнорируется. "input option alredy specificated" - спецификация INPUT уже сделана, эта игнорируется. "input-table option alredy specificated" - спецификация INTAB уже сделана, эта игнорируется. "Extern variables option already specificated" - спецификация EXTRN уже сделана. Эта игнорируется. "There must be .K or .J suffix" - В этом месте должен быть временной суффикс .K или .J "No time suffix. Set default .J" - отсутствует временной суффикс, ставится по умолчанию .J "There must be .JK or .KL suffix" - В этом месте должен быть временной суффикс .KL или .JK "No time suffix. Set default .JK" - отсутствует временной суффикс, ставится по умолчанию .JK "There must be no time suffix" - в этом месте не должно быть временного суффикса.