Осуществление заданных движений динамических систем с минимальной энергией

Министерство образования Республики Беларусь Учреждение образования

«Гомельский государственный університет имени Франциска Скорины»

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

Дипломная работа

Осуществление заданных движений динамических систем с минимальной энергией

Исполнитель студент группы ПМ — 52

Д.А. Синица Научный руководитель кандидат физико-математических наук, А. В. Лубочкин Гомель 2013

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

Цель настоящей работы — описать подход по использованию методов оптимального управления для задачи следящих систем.

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

1. Обзор литературы

Согласно [1], «Автоматическим регулированием называется процесс поддержания или изменения по заданным условиям какой-нибудь величины в машинах, аппаратах или иных технических устройствах». Построение первого типа процессов называют задачей регулирования, второго — задачей слежения.

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

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

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

2. Постановка задачи осуществления заданных движений

Пусть на промежутке времени динамическая система с управлением описывается уравнением

(2.1)

Где

п — вектор состояния системы в момент ,

значение скалярного управляющего воздействия,

.

Будем считать, что доступные управления стеснены неравенством

.

Наряду с уравнением (2.1) рассмотрим движение на фазовой плоскости

(2.2)

заданное кусочно-гладкой функцией .

Будем говорить, что движение (2.2) допустимо (осуществимо), если существует такое доступное управление

Что Пусть — область фазового пространства системы (2.1), внутренность которой содержит движение (2.2):

Определение. Функция

(2.3)

называется ограниченной обратной связью, осуществляющей движение (2.2), если

1)

2)

3) замкнутая система

(2.4)

имеет решение ;

4) решение системы (2.4) асимптотически устойчиво.

При фактическом построении обратной связи для следящих систем стараются дополнительно обеспечить следующие свойства;

5) область притяжения G достаточно большая;

6) переходный процесс обладает хорошими в определенном смысле характеристиками.

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

Следуя концепции Ляпунова из классической теории устойчивости о невозмущенном возмущенном движениях, приведем другую (эквивалентную) постановку задачи осуществления движений.

Опираясь на уравнение (2.1) (невозмущенного движения) и заданное движение, введем новые переменные Их поведение подчиняется уравнению (возмущенного движения)

(2.5)

и неравенствам ;

В результате, задача устойчивого осуществления движения (2.2) системой (2.1) с постоянным ограничением на управление свелась к задаче стабилизации тривиального решения системы (2.5) управлениями с переменными (во времени) ограничениями.

Замечания:

1) Если движение (2.2) представляет состояние равновесия с управлением то задача осуществления движения называется классической задачей регулирования и в силу приведенных операций сводится к классической задаче стабилизации с постоянными ограничениями на управление.

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

3) В задача осуществления движения рассматривается в более общей постановке, когда в классе доступных управлений не существует управлений осуществляющих движение В этом случае ставится задача приближенного осуществления движения. Если — периодическое движение, то новую задачу можно свести к рассмотренной выше. Для этого сначала в классе доступных управлений находится такое управлениеt> 0, и начальное состояние, что соответствующая им Т — периодическая траектория системы (2.1) наилучшим образом приближает движение :

(2.6)

Затем по описанной в данной работе схеме строится обратная связь, осуществляющая движение (норму в (2.6) можно взять другую).

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

Обратную связь (2.3) назовем дискретной (с периодом квантования), если порождаемая им траектория замкнутой системы (2.4) с начальным условием

строится по следующему правилу

Таким образом, при использовании дискретной обратной связи, отпадает вопрос о существовании решения замкнутой системы (2.4).

Очевидны изменения, которые нужно внести в определение осуществимого движения.

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

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

3. Сопровождающая линейно-квадратичная задача оптимального управления

Функцию назовем дискретным управлением с периодом квантования х > 0, если Пусть — произвольный момент времени. В классе дискретных управлений рассмотрим задачу оптимального управления.

(3.1)

где — текущий момент времени.

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

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

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

(3.2)

будем называть оптимальным (стартовым) управлением типа обратной связи.

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

В качестве обратной связи (2.3), решающей задачу осуществления движения, возьмем функцию (3.2):

(3.3)

Покажем, что функция обладает свойствами 1) — 4).

4. Свойства оптимальной стартовой обратной связи

Из (3.3) видно, что т. е. — осуществимое движение:

По определению дискретной обратной связи, замкнутая система (2.4) с и имеет решение при любом .

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

Функцию, , возьмем в качестве функции Ляпунова. Ясно, что, ,, , , — непрерывная функция. Докажем, что функция, , убывает вдоль каждой траектории системы Пусть — произвольное текущее состояние системы, - соответствующее ему значение функции Ляпунова.

Используя формулу Коши

исключим из задачи (3.1) переменные состояния и запишем ее с учетом используемого класса доступных управлений в эквивалентной функциональной форме:

(4.1)

где ,—фундаментальная матрица решений однородной системы

Пусть в произвольный текущий момент система (2.4) оказалась в состоянии), которое соответствует произвольному начальному состоянию

и реализации. Предположим, что для состояния

построен оптимальный план

задачи (4.1). На нем критерий качества задачи (4.1) принимает значение

В момент под действием управления

система (2.4) окажется в состоянии

Вектор

является планом задачи (4.1), если в ней вектор

заменить на

.

Значение критерия качества задачи (4.1) на этом плане равно

Если то

При имеем

и план будет оптимальным для задачи (4.1).

Поскольку

при следующем шаге в силу чего получим

Таким образом, функция

,

убывает вдоль последовательности Предельным значением может быть только 0. Используя ограниченность управлений и рассуждения, типичные для метода функций Ляпунова [5], нетрудно показать, что

Качество переходных процессов в замкнутой системе определяется следующим экстремальным свойством:

Если значение трактовать как затрату энергии в момент времени на переходной процесс из состояния, то последнее неравенство означает, что затрату энергии на весь переходной процесс не превосходит минимальной ее затраты на перевод за время системы (2.1) из состояния в состояние

5. Алгоритм построения обратной связи

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

(5.1)

критерия оптимальности с учетом формул

(5.2)

(5.3)

(5.4)

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

(5.5)

(5.6)

(5.7)

(5.5)-(5.7):

1) опорные компоненты, «отвечают» только за соблюдение основного ограничения (5.6) и могут принимать любые значения из отрезков:

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

2) неопорные компоненты, делятся на 2 группы: одни из них принимают критические значения (или) в сочетании с определенным знаком оценок, (см.(5.1)); вторая часть принимает некритические значения с одновременным равенством нулю их оценок. При этом очевидно, что поскольку оценки зависят от плана, то любой малый сдвиг от точки сразу нарушит равенство нулю оценок компонент второй группы (некритических компонент), т. е. сразу будет нарушена третья группа соотношений из (5.1). Если при этом у компонент, принимающих критические значения, оценки ненулевые, то малый сдвиг от точки не изменит знаки этих оценок, т. е. при малом сдвиге две первые группы соотношений из (5.1) будут выполняться.

Таким образом, можно выделить две группы оптимальных признаков плана задачи (5.5)-(5.7):

1) устойчивые признаки:

(5.8)

2) неустойчивые признаки:

(5.9)

На вышеизложенном анализе структуры оптимального плана задачи (5.5)-(5.7), на разделении оптимальных признаков на две группы (5.8) и (5.9) и базируются прямые методы решения этой задачи.

Направление улучшения плана в случае невыполнения соотношений оптимальности (5.1), т. е. в случае невыполнения оптимальных признаков (5.8), (5.9), подсказывает доказательство необходимой части критерия оптимальности.

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

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

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

Описываемые прямые методы (метод первого порядка и метод второго порядка) различаются множествами накапливаемых оптимальных признаков и, в соответствии с этим — способами накопления этих признаков.

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

Метод второго порядка постепенно накапливает все оптимальные признаки. Для удержания неустойчивых признаков (5.9) привлекается специальная конструкция — опора целевой функции. Таким образом, метод второго порядка (в отличие от метода первого порядка) кроме опоры ограничений использует также опору целевой функции.

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

Теперь рассмотрим сам алгоритм. Пусть известны: начальное состояние системы (2.1), движение и управление До начала процесса функционирования системы в момент при выбранных значениях, построим оптимальное программное управление, задачи (4.1), соответствующее начальному состоянию . Это можно сделать, например, методами линейного программирования [9−10], так как все элементы задачи (4.1) известны. В момент, когда начался процесс слежения, управление подается на вход системы (2.1), приводит ее в момент в состояние

Предположим, что в моменты 0,, …, вычислены значения

…,

обратной связи (3.3) и система (2.1), замкнутая обратной связью (3.3), оказалась в состоянии По предположению, в предыдущий момент когда система находилась в состоянии задача (4.1) с.

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

Двойственный метод построит оптимальную опору (и соответствующее ей оптимальное управление. Построив решение

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

(5.10)

(5.11)

(5.12)

Здесь — симметрическая неотрицательно определенная — матрица; - векторы;

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

Пусть — матрица основных ограничений задачи, опора задачи (5.10)-(5.12). Псевдопланом задачи (3.1) назовем — вектор, удовлетворяющий ограничениям

(5.13)

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

(5.14)

Компоненты псевдоплана,, назовем опорными, , — неопорными. Положим Согласованный опорный псевдоплан назовем невырожденным, если множество является опорой целевой функции (5.10). В противном случае согласованный опорный план будем называть вырожденным.

Критерий оптимальности. Для того, чтобы согласованный опорный псевдоплан был решением задачи (5.10)-(5.12), достаточно, а в случае его невырожденности и необходимо, вычисление соотношений:

(5.15)

где .

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

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

Программа DMQP — программа конечного двойственного метода линейно-квадратичного программирования.

Как и прежде рассматриваем задачу линейно-квадратичного программирования:

где — двумерный массив основных ограничений размерности ,

и — векторы-столбцы размерности (в нашем случае),

и — нижнее и верхнее ограничениям на переменные,

— вектор-столбец размерности, соответствующий начальному плану задачи,

— количество основных ограничений задачи,

— количество переменных.

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

Входной файл имеет структуру представленую в таблице 1.

Таблица 1 — Структура входного файла

В файле:

— транспонированный двумерный массив ,

и соответственно равны и ,

и соответственно равны и ,

— размер опоры ограничений начального согласованного опорного плана,

— обратная матрица к опорной матрице,

— размер опоры целевой функции начального согласованного опорного плана,

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

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

KOD (2.4) — код завершения подпрограммы, принимает значение:

KOD (2.4)=0 — полученный оптимальный план задачи,

KOD (2.4)=1 — исчерпан ресурс итераций,

KOD (2.4)=2 — получено направление бесконечного убывания целевой функции (массив R),

KOD (2.4)=3 — вычислительный процесс остановлен из-за накопления ошибок округлений (точка, содержащаяся в массиве X, не удовлетворяет ограничениям задачи),

KOD (2.4)=4 — вычислительный процесс остановлен из-за недостаточных размеров массива E (E — двумерный массив содержит при в первых строках и первых столбцах опорную матрицу), так как количество индексов, вошедших в опору ограничений, превысило значение первого измерения, двухмерного массива E.

Выходной файл содержит следующую результирующую информацию:

— значение параметра KOD (2.4). Если он равен нулю, то построен оптимальный план задачи;

— размерность опоры ограничений;

— опору ограничений;

— размерность опоры целевой функции;

— опору целевой функции;

— неопорные индексы;

— оптимальный план задачи;

6. Примеры

Пример 1. Рассмотрим систему управления

(6.1)

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

(6.2)

и построим обратную связь (3. 3), с которой замкнутая система

(6.3)

имеет периодическое движение (6.2) в качестве асимптотически устойчивого предельного цикла.

При построении обратной связи были выбраны следующие параметры вспомогательной задачи оптимального управления (4.1):

Результаты решения задачи:

Рисунок 1 и 2 приведены для движения с началом в точке (-2,0)

Рисунок 1 — Реализация управления

На рисунке видно, что фунция стабилизируется в 0.

Рисунок 2 — Изменение функции Ляпунова

Рисунок 3 - Фазовые траектории замкнутой системы.

С началом движения в точках:

— 1 — (-2,0);

— 2 — (-0.5,0);

— 3 — (0.5,0);

— 4 — (2,0);

— 5 — траектория предельного цикла.

Пример 2. Рассмотрим систему управления

(6.4)

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

(6.5)

С помощью решения вспомогательной задачи (4.1) при

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

Результаты решения задачи:

Рисунок 1 и 2 приведены для движения с началом в точке (-2,0)

Рисунок 4 — Реализация управления

На рисунке видно, что фунция стабилизируется в -2.

Рисунок 5 — Изменение функции Ляпунова

Рисунок 6 - Фазовые траектории замкнутой системы.

С началом движения в точках:

— 1 — (2,0);

— 2 — (2,1.5);

— 3 — (2,2.5);

— 4 — (2,4);

— 5 — траектория предельного цикла.

Заключение

оптимальный управление дискретный квадратичный

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

Список использованных источников

1 Айзерман, М. А. Лекции по теории автоматического регулирования. М.: Физматгиз, 1958.

2 Габасов, Р., Кириллова, Ф.М., Ружицкая, Е. А. Решение классической задачи регулирования методами оптимального управления // АиТ. 2001. М. С. 18−29.

3 Барбашип, Е.А.

Введение

в теорию устойчивости. М.: Наука, 1967.

4 Квакерпаак X., Сиван Р. Линейные оптимальные системы управления. М.: Мир, 1977.

5 Малкин И. Г. Теория устойчивости движения. М.: Наука, 1966.

6 Поптрягин, Л.С., Болтянский, В.Г., Гамкрелидзе, Р.В., Мищенко, Е. Ф. Математическая теория оптимальных процессов. М.: Наука, 1969.

7 Габасов, Р., Кириллова, Ф.М., Костюкова, О. И. Построение оптимальных управлений типа обратной связи в линейной задаче // Доклады АН СССР, 1991. Т.320. Ж.С.1294−1299.

8 Gabasov, R., Kirillova, F.M., Ruzhitskaya, Е.А. Stabilization of dynamical systems with the help of optimization methods // Nonsmooth and discontinuous problems of control and optimization. Proceedings volume from the IFAC Workshop, Chelyabinsk, Russia, 1998.P.35−41.

9 Габасов, P., Кириллова, Ф.М., Тятюшкин, А. И. Конструктивныеметодыоптимизации. Часть 1. Линейные задачи. Минск: Изд-во «Университетское», 1984.

10 Данциг, Дж. Линейное программирование, его применения и обобщения. М.: Прогресс, 1966.

11 Габасов, Р., Кириллова, Ф.М., Ружицкая, Е. А. Гашение колебаний струны // ДНАН Беларуси. 1999. Т.43, № 5, С.9−12.

12 Габасов, Р. Конструктивные методы оптимизации: в 5 ч. Ч. 1. Линейные задачи / Р. Габасов, Ф. М. Кириллова, А. И. Тятюшкин. — Мн.: БГУ, 1983. — 214 с.

13 Габасов, Р. Конструктивные методы оптимизации: в 5 ч. Ч. 2. Задачи управления / Р. Габасов, Ф. М. Кириллова. — Мн.: Университетское, 1984. — 204 с.

14 Конструктивные методы оптимизации: в 5 ч. Ч. 4. Выпуклые задачи / Р. Габасов [и др.]. — Мн.: Университетское, 1987. — 223 c.

15 Ракецкий, В. М. Решение общей задачи выпуклого квадратичного программирования прямым методом / В. М. Ракецкий // Програм. обеспечение ЭВМ. — Мн.: Ин-т математики АН БССР, 1985. — Вып. 55. — С. 113−123.

16 Ракецкий, В. М. Решение общей задачи выпуклого квадратичного программирования двойственным методом / В. М. Ракецкий // Програм. обеспечение ЭВМ. — Мн.: Ин-т математики АН БССР, 1985. — Вып. 55. — С. 124−129.

17 Габасов, Р. Построение оптимальных управлений типа обратной связи в линейной задаче / Р. Габасов, Ф. М. Кириллова, О. И. Костюкова // Докл. АН СССР. — 1991. — Т. 320,№ 6. — С. 1294−1299.

18 Габасов, Р. К методам стабилизации динамических систем / Р. Габасов, Ф. М. Кириллова, О. И. Костюкова // Изв. РАН. Техн. кибернетика. — 1994. -№ 3. — С. 67−77.

19 Габасов, Р. Стабилизация линейных динамических систем оптимальными управлениями линейно-квадратичных задач / Р. Габасов, А. В. Лубочкин // Прикл. матем. и механ. — 1998. — Т. 62, Вып. 4. — С. 556−565.

20 Лубочкин, А. В. Использование ограниченных оптимальных обратных связей линейно-квадратичных задач для осуществления заданных движений динамических систем / А. В. Лубочкин // IX Белорусская матем. конф.: тез. докл.: в 3 ч. Ч. 3 (2004;Гродно), 3−6 ноября 2004 г., г. Гродно. — Гродно: ГрГУ, 2004. — С. 124−125.

Приложения

Приложение А

Текст программы 1

#pragma hdrstop

# include

# include

# include

# include

# include

# include

#define M 2

#define N 25

#define N_M 23

#define TBEG0.0

#define TETA1.0

#define L 4.0

#define NX 2

double tau,

x[NX] = { 0.5, 0.0 },

z0[NX] = { 1.0, 0.0 },

z[NX],// Точка на пред. цикле y=y (t):

// z (t) = y (t+TETA) = F (t+TETA)*z0

b[M], A[M][N], dn[N], dw[N], eps1 = 0.1, eps2=0.1,

P[M][M], **Pr, Q[M][M], **D, v[N_M], **Ma, u[N], Del[N], J, a1, a2;

int Jop[M], Jn[N_M], kn,

Joc[N_M], koc, Jnn[N_M], knn,

kod4;

void callocPr (int);

void freePr (int);

void puN1NIntU2 (double);

void puIntU2 (double);

void formParam (double);

void formElem (void);

void procsogl (void);

int gs (int, double**, /* Метод Гаусса */

double*, int*);

void f_x (double, double*,/* f1(t) = x2(t) */

double*); /* f2(t) = -x1(t) + u (t) */

void FEHL (int, double,

double*, void (f_p)(double, double*, double*),

double*, double*,

double*, double*,

double*, double*,

double*, double*);

void RKF45 (int, double, /* Метод Рунге-Кутта-Фельдберга */

double, double*,

void (f_p)(double, double*, double*), double*,

double*, double*,

double*, double*,

double*, double*,

int*, int*,

int*, int*,

int*, int*,

double*, double*,

double*, double*);

void prifl0 (int*, double*, double*);

FILE *d;

FILE *r;

FILE *ru; /* Файл. dat: t, u (t) */

FILE *rJ; /* Файл. dat: t, J (t) */

FILE *rx; /* Файл. dat: t, x (t) */

#pragma argsused

int main (int argc, char* argv[])

{

char f_Name_u[13],/* Имя файла. dat: t, u (t) */

f_Name_J[13],/* Имя файла. dat: t, J (t) */

f_Name_x[13]; /* Имя файла. dat: t, x (t) */

double h, t, tout, ypx [NX], hx=0, f1x [NX], f2x [NX], f3x [NX], f4x [NX], f5x [NX], rerr = 1.0e-12, aerr = 0.0,

savrex = 0, savaex = 0;

int j, Ostanov, inte, iflx, jflx, kflx, nfex, kopx, initx;

D = (double**)calloc (N_M, sizeof (double));

for (j=0; j

D[j] = (double*)calloc (N_M, sizeof (double));

Ma = (double**)calloc (N_M, sizeof (double));

for (j=0; j

Ma[j] = (double*)calloc (N_M, sizeof (double));

/* Задание возмущения: */

printf («Возмущение: a1 * sin (a2 * t) n»);

printf («a1=»); scanf («%lf», &a1);

printf («a2=»); scanf («%lf», &a2);

/* Открытие файлов: */

printf («nFileName: »); scanf («%s», f_Name_u);

strcpy (f_Name_x, f_Name_u);

strcpy (f_Name_J, f_Name_u);

strcat (f_Name_u, «_u.dat»);

strcat (f_Name_J, «_J.dat»);

strcat (f_Name_x, «_x.dat»);

if ((ru = fopen (f_Name_u, «wt»)) == NULL)

{

perror (f_Name_u);

exit (-1);

}

if ((rJ = fopen (f_Name_J, «wt»)) == NULL)

{

perror (f_Name_J);

exit (-1);

}

if ((rx = fopen (f_Name_x, «wt»)) == NULL)

{

perror (f_Name_x);

exit (-1);

}

tau = TBEG;

/* Печать в файл и на экран x0: */

fprintf (rx, «%lf «, tau);

printf («ntau = %lf «, tau);

printf («x = («);

for (j=0; j

{

fprintf (rx, «%lf «, x[j]);

printf («%lf «, x[j]);

}

fprintf (rx, «n»);

printf («)n»);

/* Начальн. формирование некотор. параметров конечномерн.задачи: */

h = TETA/N; /* шаг (длина перем.)*/

printf («h = %lfn», h);

getchar ();

/* Работа регулятора: */

Ostanov = 0;

while (!Ostanov)

{

/* Программн. упpавл. задачи с перем. структурой: */

puN1NIntU2 (h);

/* Печ. в файл u (t), x (t) на [tau, tau+h]: */

fprintf (ru, «%lf %lfn», tau, u[0]);

fprintf (ru, «%lf %lfn», tau+h, u[0]);

fprintf (rJ, «%lf %lfn», tau, J);

tout = tau;

printf («tau = %lfn», tau);

while (tout < tau+h-0.1)

{

t = tout;

if (tau + h — tout > 0.1)

tout += 0.1;

else

tout = tau + h;

printf («nt = %lf», t);

printf («tout = %lf», tout);

iflx = 1;

inte = 1;

while (inte)

{

RKF45 (NX, t, tout, x, f_x, &hx, ypx, f1x, f2x, f3x, f4x, f5x,

&iflx, &nfex, &kopx, &initx, &jflx, &kflx,

&rerr, &aerr, &savrex, &savaex);

if (iflx==2)

inte = 0;

else

{

printf («xn»);

prifl0 (&iflx, &rerr, &aerr);

}

}

fprintf (rx, «%lf «, tout);

printf («x = («);

for (j=0; j

{

fprintf (rx, «%lf «, x[j]);

printf («%lf «, x[j]);

}

fprintf (rx, «n»);

printf («) »);

}

tau += h;

printf («nnn tau = %lfnnn», tau);

while (1)

(Ostanov==1))

break;

} /* while (!Ostanov) */

fcloseall ();

for (j=0; j

free (D[j]);

free (D);

for (j=0; j

free (Ma[j]);

free (Ma);

return 0;

}

void callocPr (int n)

{

int i;

Pr = (double**)calloc (n, sizeof (double));

for (i=0; i

Pr[i] = (double*)calloc (n, sizeof (double));

}

void freePr (int n)

{

int i;

for (i=0; i

free (Pr[i]);

free (Pr);

}

void puN1NIntU2 (double h)/* Прогр.упр.задачи пер. стр.*/

{

int i;

/* Программн. упр-ие миним. энергии: */

kod4 = 0;

puIntU2 (h);

if (kod4==0)

{

J = 0.0;

for (i=0; i

J += u[i]*u[i]/2.0;

J *= h;

printf («J = %.12lfn», J);

}

else

{

printf («Нет решения!!! n»);

getchar ();

}

} /* puN1NIntU2 */

void puIntU2 (double h)/* Прогр.упр.мин.эн. */

{

double ur,

**H;

int i, j, k,

kod,

nn,

UslSogl, UslOpt;

H = (double**)calloc (N_M, sizeof (double));

for (j=0; j

H[j] = (double*)calloc (N_M, sizeof (double));

/* Формирование параметров конечномерн. задачи: */

formParam (h);

/* Формирование элементов опорного псевдоплана: */

formElem ();

/* Пров-ка согласов-ти оп. пс-плана: */

UslSogl = 1;

for (j=0; j

if (((Del[ Jnn[j] ] < -eps1) &&

(fabs (u[ Jnn[j] ] - dn[ Jnn[j] ]) < eps2)) ||

((Del[ Jnn[j] ] > eps1) &&

(fabs (u[ Jnn[j] ] - dw[ Jnn[j] ]) < eps2)))

{

UslSogl = 0;

break;

}

printf («UslSogl = %dn», UslSogl);

if (!UslSogl)

procsogl ();

UslOpt = 1;

for (j=0; j

if ((u[Jop[j]] < dn[Jop[j]] - eps2) ||

(u[Jop[j]] > dw[Jop[j]] + eps2))

{

UslOpt = 0;

break;

}

if (!UslOpt)

{

/* Реш.конечномерн.з-чи дв. методом: */

/* Открытие файла для передачи в прогр. *.for (двойств.метод) */

if ((d = fopen («dat», «wt»)) == NULL)

{

perror («dat»);

exit (-1);

}

/* Выч.обр.м. к м. D */

callocPr (koc);

for (j=0; j

{

for (i=0; i

{

for (k=0; k

Pr[i][k] = D[i][k];

v[i] = 0.0;

}

v[j] = 1.0;

gs (koc, Pr, v, &kod);

if (kod)

{

printf («nОпорная матрица (цел.ф-ии) вырождена!!! n»);

getchar ();

}

for (i=0; i

H[i][j] = v[i];

}

freePr (koc);

/* Печать параметров задачи в файл */

fprintf (d, «%d %dn», M, N);/* размеры M, N */

for (j=0; j

{

for (i=0; i

fprintf (d, «%lf «, A[i][j]);

fprintf (d, «n»);

}

for (i=0; i

fprintf (d, «%lf «, b[i]);

fprintf (d, «n»);

for (i=0; i

fprintf (d, «%lf %lfn», dn[i], dw[i]);

fprintf (d, «%dn», M);/* размер оп.огр. M */

for (i=0; i

{

for (j=0; j

fprintf (d, «%lf «, Q[i][j]);

fprintf (d, «n»);

}

fprintf (d, «%dn», koc); /* разм.оп.ц.ф. koc */

for (j=0; j

for (i=0; i<=j; i++)

fprintf (d, «%lfn», H[i][j]);

for (i=0; i

fprintf (d, «%dn», Jop[i] + 1);

for (i=0; i

fprintf (d, «%dn», Joc[i] + 1);

for (i=0; i

fprintf (d, «%dn», Jnn[i] + 1);

for (j=0; j

fprintf (d, «%lfn», u[j]);

fclose (d);

system («errfile.exe»);

system («dmqp.exe»);

/* Открытие файла для чтения из прогр. *.for (двойств.метода) */

if ((r = fopen («res», «r»)) == NULL)

{

perror («res»);

exit (-1);

}

fscanf (r, «%d», &kod4);

printf («Результаты:nkod (4) = %dn», kod4);

if (kod4==0)

{

printf («Jop = («);

fscanf (r, «%d», &i);

nn = i;

for (j=0; j

{

fscanf (r, «%d», &k);

Jop [j] = k — 1;

printf («%d «, Jop[j]);

}

printf («)n»);

printf («Joc = («);

fscanf (r, «%d», &i);

nn += i;

for (j=0; j

{

fscanf (r, «%d», &k);

Joc[j] = k — 1;

printf («%d «, Joc[j]);

}

printf («)n»);

printf («Jnn = («);

for (j=0; j

{

fscanf (r, «%d», &k);

Jnn[j] = k — 1;

printf («%d «, Jnn[j]);

}

printf («)n»);

for (k=0; k

{

fscanf (r, «%lf», &ur);

u[k] = ur;

}

printf («u =n»);

for (i=0; i<10; i++)

printf («%7.1lf», u[i]);

printf («n»);

for (i=10; i<20; i++)

printf («%7.1lf», u[i]);

printf («n»);

for (i=20; i

printf («%7.1lf», u[i]);

printf («n»);

} /* if (kod4==0) */

fclose ®;

} /* if (!UslOpt) */

for (j=0; j

free (H[j]);

free (H);

} /* puIntU2 */

void formParam (double h)/* Форм.парам.конечном.з-чи */

{

double t_tn, t_tw,

te, tauNh;

int i, j;

for (j=0; j

{

t_tn = N*h — j*h;

t_tw = N*h — (j+1)*h;

A[0][j] = cos (t_tw) — cos (t_tn);

A[1][j] = sin (t_tn) — sin (t_tw);

}

te = N*h;

tauNh = tau + N*h;

z[0] = cos (tauNh)*z0[0] + sin (tauNh)*z0[1];

z[1] = -sin (tauNh)*z0[0] + cos (tauNh)*z0[1];

b[0] = z[0] - cos (te)*x[0] - sin (te)*x[1];

b[1] = z[1] + sin (te)*x[0] - cos (te)*x[1];

for (i=0; i

{

dn[i] = -L;

dw[i] = L;

}

} /* formParam */

void formElem (void) /* Форм. эл-тов оп. пс-плана */

{

int i, j, k,

jr,

kod,

PrOp;

/* Опора огр. Jop */

Jop[0] = 0;

Jop[1] = N-1;

kn = N — M; /* неопорн.инд. Jn */

k = 0;

for (j=0; j

{

PrOp = 0;

for (i=0; i

if (j==Jop[i])

{

PrOp = 1;

break;

}

if (!PrOp)

{

Jn[k] = j;

k++;

}

}

koc = 0;

knn = kn; /* дв.неоп.инд. Jnn */

for (j=0; j

Jnn[j] = Jn[j];

for (i=0; i

for (j=0; j

P[i][j] = A[i][Jop[j]];

callocPr (M);

for (j=0; j

{

for (i=0; i

{

for (k=0; k

Pr[i][k] = P[i][k];

v[i] = 0.0;

}

v[j] = 1.0;

gs (M, Pr, v, &kod);

if (kod)

{

printf («nОпорная матрица (ограничений) вырождена!!! n»);

getchar ();

}

for (i=0; i

Q[i][j] = v[i];

}

freePr (M);

/* псевдоплан u */

for (j=0; j

u[Jnn[j]] = dn[Jnn[j]];

if (koc>0) /* uoc */

{

/* uoc */

}

if (koc>0) /* U=(*, Uoc, Unn) */

for (i=0; i

u[Joc[i]] = v[i];

callocPr (M);

for (i=0; i

{

for (j=0; j

Pr[i][j] = P[i][j];

v[i] = b[i];

for (j=0; j

v[i] -= A[i][Jn[j]] * u[Jn[j]];

}

gs (M, Pr, v, &kod);

freePr (M);

if (kod)

{

printf («nОпорная матрица (ограничений) вырождена!!! n»);

getchar ();

}

for (i=0; i

u[Jop[i]] = v[i];

callocPr (M);

for (i=0; i

{

for (k=0; k

Pr[i][k] = P[k][i];

v[i] = u[Jop[i]];

}

gs (M, Pr, v, &kod);

freePr (M);

if (kod)

{

printf («nОпорная матрица (ограничений) вырождена!!! n»);

getchar ();

}

for (j=0; j

{

Del[j] = u[j];

for (i=0; i

Del[j] -= v[i] * A[i][j];

}

} /* formElem */

void procsogl (void)

{

double l[N],

Qb[M],

Th, Thj;

int i, j, k,

kod,

iterS,

prTh,

j0;

iterS = 0;

while (1)

{

printf («iterS = %dn», iterS + 1);

/* l */

for (j=0; j

if (Del[Jnn[j]] > eps1)

l[Jnn[j]] = dn[Jnn[j]] - u[Jnn[j]];

else

if (Del[Jnn[j]] < -eps1)

l[Jnn[j]] = dw[Jnn[j]] - u[Jnn[j]];

else

l[Jnn[j]] = 0.0;

if (koc>0) /* loc */

{

callocPr (M);

for (j=0; j

{

for (i=0; i

{

for (k=0; k

Pr[i][k] = P[i][k];

v[i] = A[i][Joc[j]];

}

gs (M, Pr, v, &kod);

if (kod)

{

printf («nОпорная матрица (ограничений) вырождена!!! n»);

getchar ();

}

for (i=0; i

Ma[i][j] = v[i];

}

freePr (M);

for (i=0; i

for (j=0; j

{

D[i][j] = 0.0;

for (k=0; k

D[i][j] += Ma[k][i] * Ma[k][j];

}

for (i=0; i

D[i][i] += 1.0;

callocPr (M);

for (i=0; i

{

for (k=0; k

Pr[i][k] = P[i][k];

Qb[i] = 0.0;

for (j=0; j

Qb[i] += A[i][Jnn[j]] * l[Jnn[j]];

}

gs (M, Pr, Qb, &kod);

freePr (M);

if (kod)

{

printf («nОпорная матрица (ограничений) вырождена!!! n»);

getchar ();

}

for (i=0; i

{

v[i] = 0.0;

for (j=0; j

v[i] -= Ma[j][i] * Qb[j];

}

callocPr (koc);

for (i=0; i

for (k=0; k

Pr[i][k] = D[i][k];

gs (koc, Pr, v, &kod);

freePr (koc);

if (kod)

{

printf («nОпорная матрица (цел.ф-ии) вырождена!!! n»);

getchar ();

}

}

for (j=0; j

l[Joc[j]] = v[j];

callocPr (M);

for (i=0; i

{

for (j=0; j

Pr[i][j] = P[i][j];

v[i] = 0.0;

for (j=0; j

v[i] -= A[i][Jn[j]] * l[Jn[j]];

}

gs (M, Pr, v, &kod);

freePr (M);

if (kod)

{

printf («nОпорная матрица (ограничений) вырождена!!! n»);

getchar ();

}

for (j=0; j

l[Jop[j]] = v[j];

/* tnn */

callocPr (M);

for (i=0; i

{

for (k=0; k

Pr[i][k] = P[k][i];

Qb[i] = l[Jop[i]];

}

gs (M, Pr, Qb, &kod);

freePr (M);

if (kod)

{

printf («nОпорная матрица (ограничений) вырождена!!! n»);

getchar ();

}

for (j=0; j

{

v[j] = l[Jnn[j]];

for (i=0; i

v[j] -= Qb[i] * A[i][Jnn[j]];

}

Th = 1; /* Theta */

prTh = -1;

for (j=0; j

if (l[Joc[j]] < -eps2)

{

Thj = (dn[Joc[j]] - u[Joc[j]]) / l[Joc[j]];

if (Thj < Th)

{

Th = Thj;

j0 = j; /* ??? */

prTh = 0;

}

}

else

if (l[Joc[j]] > eps2)

{

Thj = (dw[Joc[j]] - u[Joc[j]]) / l[Joc[j]];

if (Thj < Th)

{

Th = Thj;

j0 = j;

prTh = 0;

}

}

for (j=0; j

if (Del[Jnn[j]] * v[j] < -eps1)

{

Thj = -Del[Jnn[j]] / v[j];

if (Thj < Th)

{

Th = Thj;

j0 = j;

prTh = 1;

}

}

for (j=0; j

u[j] += Th*l[j];

if (Th == 1)

{

break;

}

else

if (prTh)

{

Joc[koc] = Jnn[j0];

koc++;

knn—;

for (i=j0; i

Jnn[i] = Jnn[i+1];

}

else

{

Jnn[knn] = Joc[j0];

knn++;

koc—;

for (i=j0; i

Joc[i] = Joc[i+1];

}

callocPr (M);

for (i=0; i

{

for (k=0; k

Pr[i][k] = P[k][i];

v[i] = u[Jop[i]];

}

gs (M, Pr, v, &kod);

freePr (M);

if (kod)

{

printf («nОпорная матрица (ограничений) вырождена!!! n»);

getchar ();

}

for (j=0; j

{

Del[j] = u[j];

for (i=0; i

Del[j] -= v[i] * A[i][j];

}

iterS++;

} /* while (1) */

}

int gs (int m, double **pX,

double copX[], int *kod)

{

double tol = 1E-5,

bigpX, pXpX,

sav;

int i, j, k, ky, imax,

ix, jx, i0, ib;

for (j=0; j

{

bigpX = 0.;

for (i=j; i

{

pXpX = fabs (pX [i] [j]);

if (fabs (bigpX) < pXpX)

{

bigpX = pX [i] [j];

imax = i;

}

}

if (fabs (bigpX) <= tol)

{

*kod = 1;

return 1;

}

for (k=j; k

{

sav = pX [imax] [k];

pX [imax] [k] = pX [j] [k];

pX [j] [k] = sav / bigpX;

}

sav = copX [imax];

copX [imax] = copX [j];

copX [j] = sav / bigpX;

if (j==m-1)

break;

for (ix=j+1; ix < m; ix++)

{

for (jx=j+1; jx < m; jx++)

pX [ix] [jx] -= pX [ix] [j] * pX [j] [jx];

copX [ix] -= copX [j] * pX [ix] [j];

}

}

for (ky=0; ky < m-1; ky++)

{

ib = m-2-ky;

i0 = m-1;

for (k=0; k <= ky; k++)

{

copX [ib] -= pX [ib] [i0] * copX [i0];

i0—;

}

}

*kod=0;

return 0;

}

void f_x (double t, /* Выч. f1(t) = x2(t) */

double y[], /* f2(t) = -x1(t) + u (t) */

double f[])

{

f[0] = y[1];

f[1] = -y[0] + u[0] + a1*sin (a2*t);

} /* f_x () */

# define True 1

# define False 0

void FEHL (int neqn, double t,

double y[], void (f_p)(double, double*, double*),

double *hc, double yp[],

double f1[], double f2[],

double f3[], double f4[],

double f5[], double s[])

{ double h;

double ch, rab;

int k;

void (*fun)(double, double*, double*);

h=*hc;

fun = f_p;

ch = h/4.0;

for (k=0; k

f5[k] = y[k] + ch * yp[k];

rab = t + ch;

fun (rab, f5, f1);

ch = 3.0 * h / 32.0;

for (k=0; k

f5[k] = y[k] + ch * (yp[k] + 3.0 * f1[k]);

rab = t + 3.0 * h / 8.0;

fun (rab, f5, f2);

ch = h / 2197.0;

for (k=0; k

f5[k] = y[k] + ch * (1932.0 * yp[k] + (7296.0 * f2[k] - 7200.0 * f1[k]));

rab = t + 12.0 * h / 13.0;

fun (rab, f5, f3);

ch = h / 4104.0;

for (k=0; k

f5[k] = y[k] + ch * ((8341.0 * yp[k] - 845.0 * f3[k]) +

(29 440.0 * f2[k] - 32 832.0 * f1[k]));

rab = t + h;

fun (rab, f5, f4);

ch = h / 20 520.0;

for (k=0; k

f1[k] = y[k] + ch * ((-6080.0 * yp[k] + (9295.0 * f3[k] ;

5643.0 * f4[k])) + (41 040.0 * f1[k] - 28 352.0 * f2[k]));

rab = t + h / 2.0;

fun (rab, f1, f5);

/*ВЫЧИСЛИТЬ ПРИБЛИЖЕННОЕ РЕШЕНИЕ В ТОЧКЕ t + h*/

ch = h / 7 618 050.0;

for (k=0; k

s[k] = y[k] + ch * ((902 880.0 * yp[k] + (3 855 735.0 * f3[k] ;

1 371 249.0 * f4[k])) + (3 953 664.0 * f2[k] + 277 020.0 * f5[k]));

*hc=h;

return;

}

void RKF45 (int neqn, double t,

double tout, double y[],

void (f_p)(double, double*, double*), double *hc,

double yp[], double f1[],

double f2[], double f3[],

double f4[], double f5[],

int *iflagc, int *nfec,

int *kopc, int *initc,

int *jflagc, int *kflagc,

double *relerrc, double *abserrc,

double *savrec, double *savaec)

{

int iflag, nfe, kop, init, jflag, kflag;

double h, savre, savae,

relerr, abserr;

int maxnfe = 30 000;

double remin = 1.0e-15;

int hfaild, output;

double a, ae, dt, ee, eeoet, esttol, et, hmin, eps, u26,

rer, s, scale, tol, toln, ypk, rab1;

int k, mflag;

void (*fun)(double, double*, double*);

fun = f_p;

iflag=*iflagc; nfe=*nfec; kop=*kopc; init=*initc; jflag=*jflagc;

kflag=*kflagc;

h=*hc; savre=*savrec; savae=*savaec;

relerr = *relerrc;

abserr = *abserrc;

/*Проверить входные параметры*/

if (neqn < 1)

goto m10;

if (relerr < 0.0 || abserr < 0.0)

goto m10;

mflag=abs (iflag);

if (mflag == 0 || mflag > 8)

goto m10;

if (mflag ≠ 1)

goto m20;

/*Первый вызов, вычислить машинное эпсилон*/

eps = 1.0;

while (eps+1.0 > 1.0)

eps = eps/2.0;

u26 = 26.0*eps;

goto m50;

/*Ошибки во входной информации*/

m10: iflag = 8;

goto mexit;

/*Проверить возможность продолжения*/

m20: if (t==tout && kflag ≠ 3)

goto m10;

if (mflag ≠ 2)

goto m25;

if (kflag == 3 || init==0)

goto m45;

if (kflag == 4)

goto m40;

if (kflag == 5 && abserr == 0.0)

goto m30;

if (kflag == 6 && relerr <= savre && abserr <= savae)

goto m30;

goto m50;

m25: if (iflag == 3)

goto m45;

if (iflag == 4)

goto m40;

if (iflag == 5 && abserr > 0.0)

goto m45;

m30: goto mexit;

m40: nfe = 0;

if (mflag == 2)

goto m50;

m45: iflag = jflag;

if (kflag == 3)

mflag = abs (iflag);

m50: jflag = iflag;

kflag = 0;

savre = relerr;

savae = abserr;

rer = 2.0*eps+remin;

if (relerr >= rer)

goto m55;

/*Заданная граница относительной погрешности слишком мала*/

relerr = rer;

iflag = 3;

kflag = 3;

goto mexit;

m55: dt = tout-t;

if (mflag==1)

goto m60;

if (init==0)

goto m65;

goto m80;

m60: init = 0;

kop = 0;

a = t;

fun (a, y, yp);

nfe = 1;

if (t ≠ tout)

goto m65;

iflag = 2;

goto mexit;

m65: init = 1;

h = fabs (dt);

toln = 0.0;

for (k=0; k

{ tol = relerr*fabs (y[k])+abserr;

if (tol > 0.0)

{ toln=tol;

ypk=fabs (yp[k]);

rab1=h*h*h*h*h;

if (ypk*rab1 > tol)

h=exp (0.2*log (tol/ypk));

}

}

if (toln <= 0.0)

h=0.0;

rab1=u26*fabs (t);

if (fabs (t) < fabs (dt))

rab1=u26*fabs (dt);

if (h < rab1)

h=rab1;

jflag=-2;

if (iflag > 0)

jflag=2;

m80: if (dt < 0.0)

h=-h;

if (fabs (h) >= 2.0*fabs (dt))

kop=kop+1;

if (kop ≠ 100)

goto m85;

/*Много выходов*/

kop=0;

iflag=7;

goto mexit;

m85: if (fabs (dt) > u26*fabs (t))

goto m95;

/*Близко к выходной точке, проэкстраполировать и вернуться по месту вызова*/

for (k=0; k

y[k] = y[k] + dt * yp[k];

a=tout;

fun (a, y, yp);

nfe=nfe+1;

goto m300;

m95: output=False;

scale=2.0/relerr;

ae=scale*abserr;

/*Пошаговое интегрирование*/

m100: hfaild=False;

hmin=u26*fabs (t);

dt=tout-t;

if (fabs (dt) >= 2.0*fabs (h))

goto m200;

if (fabs (dt) > fabs (h))