Стабилизация нелинейной модели маятника с двумя нелинейностями управлениями минимальной интенсивности

 

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ

Учреждение образования

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

имени Франциска Скорины»

Математический факультет

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


Допущена к защите

Зав. кафедрой__________ А.В. Лубочкин

"____"________________20___г.



Стабилизация нелинейной модели маятника с двумя нелинейностями управлениями минимальной интенсивности

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



Исполнитель

студентка группы ПМ-52 Е. А. Полховская

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

к. ф. - м. н., доцент А. В. Лубочкин

Рецензент

доцент кафедры прикладной

математики БелГУТа,

к.т.н., доцент С. И. Жогаль





Гомель 2013

Реферат


Дипломная работа 62 страницы, 26 источников, 1 рисунок, 1 приложение.

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

Объект исследования: проблема стабилизации нелинейной модели маятника с двумя нелинейностями управлениями минимальной интенсивности.

Цель дипломной работы: метод стабилизации нелинейной модели маятника с двумя нелинейностями управлениями минимальной интенсивности.

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

Содержание


Введение

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

2. Постановка нелинейной задачи

3. Кусочно-линейная и кусочно-постоянная аппроксимация нелинейностей

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

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

. Оптимальный регулятор

. Стабилизация перевёрнутого маятника

Заключение

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

Приложение А


Введение


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

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

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

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

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


Стабилизация - одна из центральных проблем теории управления движением [1]. Классический подход к ее решению состоит в использовании разнообразных достаточных условий асимптотической устойчивости движения [2]. При этом знание подходящей функции Ляпунова позволяет оценить и качество переходных процессов [3]. Задача стабилизации является одной из причин возникновения современной теории управления и в частности теории оптимального управления. Теория оптимального управления за более чем полвека прошла огромный путь. Естественен вопрос, что дают методы оптимального управления для решения классических задач теории регулирования, актуальность которых несомненна и в наши дни. Длительное время применение теории оптимального управления к решению задач регулирования и стабилизации сдерживалось отсутствием эффективных методов синтеза оптимальных систем. Задача регулирования в классической постановке состоит в построении обратной связи, которая переводит динамическую систему из окрестности одного состояния равновесия в окрестность другого и стабилизирует систему относительно нового состояния равновесия. Первая часть задачи регулирования явилась источником возникновения современной теории оптимального управления. Вторая часть задачи исследуется в теории управления в рамках теории стабилизации. С созданием математической теории оптимальных процессов [4] появилась возможность построения оптимальных стабилизирующих обратных связей. Первый результат в этой области был получен А.М.Летовым [5] и Р.Калманом [6], которые показали, что оптимальные управления типа обратной связи для линейно-квадратичной задачи оптимального управления с бесконечным горизонтом являются линейными стабилизирующими обратными связями. Линейно-квадратичные задачи с конечным горизонтом для создания стабилизирующих обратных связей были использованы в [7-9]. Указанные примеры синтезирования стабилизирующих обратных связей на базе позиционных решений задач оптимального управления существенным образом опирались на возможность аналитического решения соответствующих задач оптимального управления и поэтому исключали прямые (геометрические) ограничения на стабилизирующие воздействия. Между тем общепризнано, что такие ограничения представляют наиболее типичный элемент современных систем управления.

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

Подход к созданию ограниченных стабилизирующих обратных связей был предложен в [10]. Он базировался на методе реализации оптимальных обратных связей, описанном в [11]. Впоследствии метод стабилизации линейных систем на базе линейных задач оптимального управления был обобщен в направлении использования позиционных решений линейно-квадратичных задач оптимального управления [12]. В дальнейшем возникла необходимость исследования задач стабилизации нелинейных систем [13, 14, 15]. В данной дипломной работе рассматривается стабилизация нелинейной модели маятника с двумя нелинейностями, которая описывается с помощью кусочно-линейной задачи минимальной интенсивности.


2. Постановка нелинейной задачи


Математический маятник - классический объект нелинейной механики. Задача его стабилизации в малом инерционным моментом на базе линейной модели и линейно-квадратичной задачи оптимального управления решена в [1]. При стабилизации в большом с помощью приложенного к оси подвеса момента линеаризованная модель уже не достаточна и приходится рассматривать нелинейное уравнение [16]


. (2.1)


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


, (2.2)


где - ускорение точки подвеса.

Уравнение (2.2) в отличие от (2.1) содержит две нелинейности и при больших возмущениях () неуправляемо.

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

В данной работе задача стабилизации динамической системы (2.2) как единая решается новым методом.

Рассмотрим динамическую систему, поведение которой в области описывается уравнением


, (2.3)


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

Функцию


(2.4)


назовем дискретной (с периодом квантования ) стабилизирующей обратной связью, если:


1);


2)траектория замкнутой системы



представляет непрерывное решение уравнения

с управлением ;


) нулевое решение , , уравнения асимптотически устойчиво и - область притяжения состояния равновесия .

Задача построения ограниченной стабилизирующей обратной связи состоит в поиске такой стабилизирующей обратной связи (2.4), которая при заданном , , удовлетворяет ограничению


.


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

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


. Кусочно-линейная и кусочно-постоянная аппроксимация нелинейностей


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

Первая процедура представляет решение задачи стабилизации для кусочно-линейной аппроксимации системы (2.3):


, (3.1)


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

Функция :


)при ;

)при проходит через точки и :

, ,

, ,

;

)при проходит через точки и (см. п.4) :

, ,

, ;

)при

(для п.3):


А - коэффициенты соответствующих прямых ломаной .

Точность аппроксимации характеризуется числом


.


Быстрые (прямые и двойственные) алгоритмы решения нескольких типов сопровождающих задач, разработанные в [19] для линейных систем, можно обобщить на кусочно-линейную систему (3.1). Это позволяет, следуя [20], реализовать оптимальные обратные связи в режиме реального времени, т.е. решить задачу стабилизации динамической системы (3.1).

Вторая процедура при решении задачи стабилизации нелинейной системы состоит в решении задачи (3.1) для кусочно-квазилинейной системы

, (3.2)


Асимптотическое решение сопровождающей задачи (5.1) для (3.2) можно строить в режиме реального времени по схеме [19], что позволит увеличить точность решения кусочно-линейной задачи (5.1), (3.1) без увеличения количества подобластей.

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

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

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


(3.3)


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

Уравнение (3.3) представляет частный случай системы


(2.1):

В качестве X возьмем область:


(3.4)


Область (3.4) разобьем, на подобласти:



Выберем следующие аппроксимации функций и :



Аппроксимации , функций представлены на рисунке 3.1 (а, б) соответственно.

a)


б)

Рисунок 3.1 - Аппроксимация нелинейностей


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


В качестве вспомогательной (сопровождающей) задачи оптимального управления в общем виде приведём задачу:


( 4.1)


где к - конкретное целое число.

А теперь рассмотрим конкретный пример.

Пусть заданы: начальное состояние , значения , , . Рассмотрим случай, когда оптимальная траектория задачи (4.1) на фазовой плоскости проходит по четырем областям: I: , II: , III: , IV: . Таким образом, маятник находится в нижнем устойчивом состоянии равновесия . Требуется поднять его в верхнее неустойчивое положение равновесия и удержать его (стабилизировать) в новом положении.

В качестве сопровождающей кусочно-линейной задачи оптимального управления (4.1) возьмем следующую задачу оптимального управления в классе дискретных управлений:


(4.2)


где


-


оптимальные моменты перехода траектории , из областей: I в II; II в III; III в IV.


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


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

Выберем число (параметр метода) и критерий качества



где если и класс доступных дискретных управлений .

Предположим, что для каждого начального состояния существует дискретное управление которое переводит систему (2.3), , в момент в состояние равновесия .

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


( 5.1)


Пусть , оптимальное программное управление задачи (5.1), - множество для которых существует решение задачи (5.1).

Функцию


(5.2)


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

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

- соответствующее оптимальное программное управление, - значение критерия качества. По определению (5.2), Траектория замкнутой системы

(5.3)

в момент окажется в состоянии , которое получается из решения уравнения



Для состояния управление ... является допустимым в задаче (5.1) и поэтому если ).

Продолжая движение вдоль траектории уравнения (5.3), найдем такое число , что



Таким образом, вдоль любой траектории системы (5.3) положительно-определенная функция , убывает. Ясно, что при . Отсюда, как при доказательстве теоремы Ляпунова об асимптотической устойчивости [2], получаем, что нулевое решение уравнения (5.3) асимптотически устойчиво в

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

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


. Оптимальный регулятор


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

Задача (4.2) задача оптимального программного управления, в ней все элементы известны, и поэтому для начального состояния она может быть решена заранее любым методом, не обращая внимание на затраты машинного времени. Однако, с точки зрения последующего синтеза, задачу (4.2) целесообразно решить адаптивным (прямым или двойственным) методом [21, 22], который кроме оптимального программного управления доставляет оптимальную опору .

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

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

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

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


. Стабилизация перевёрнутого маятника


Для стабилизации перевёрнутого маятника нам потребуется эквивалентная задача линейного программирования.

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

В первом области уравнение имеет вид:. Используя обозначения , однородная часть уравнения примет вид:


, , .


Получим фундаментальную матрицу:


, , .


Во второй области уравнение имеет вид:


.


Вводя обозначения, однородная часть уравнения примет вид:


, , .


Вычислим компоненты фундаментальной матрицы:



Получим фундаментальную матрицу:


, , , .


В третьей области уравнение имеет вид:



Если обозначим , то однородная часть уравнения примет вид:

, , .


Вычисляя



получим фундаментальную матрицу:


, , , .

В четвёртой области однородная часть уравнения имеет вид:


, ,


Фундаментальная матрица:


, , .


Пусть имеем , , , .

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


Пусть .


Используя формулу Коши и учитывая класс выбранных состояний, исключим из задачи (4.2) переменные состояния:



.


Эти значения получены путём преобразования:



Для краткости записи сделаем замену:.


Получим следующее:



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



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



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



Эти значения получены путём преобразования: для краткости записи сделаем замену:



Получаем:



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



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



Однородная часть уравнения имеет вид:



Эти значения получены путём преобразования:

Для краткости записи сделаем замену:


.


Получим следующее:



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


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



Однородная часть уравнения имеет вид:


.


Эти значения получены путём преобразования:

Для краткости записи сделаем замену:


.


Получим следующее:



Теперь в случае, когда оптимальная траектория задачи на фазовой плоскости находится в четвёртой области, получаем:

перевернутый маятник оптимальный стабилизация

Заключение


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


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


1Красовский Н.Н. Проблемы стабилизации управляемых движений / Малкин И.Г. Теория устойчивости движения. - М.: Наука, 1966.

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

3Чапаев Н.Г. О выборе параметров устойчивой механической системы // Прикладная математика и механика. - 1951. - Т.15. - Вып.2.

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

5Летов A.M. Аналитическое конструирование регуляторов // Автоматика и телемеханика. - 1960. Т.21.Ш. - С.436-441; 5. - С. 561-568, Ж. - С. 661 - 665

6Калман Р.Е. Об общей теории систем управления. - М.: АН СССР, Тр. I Конгресса ИФАК. 1961. - Т.1. - С. 521- 547.

7Kwon W.N., Pearson А.Е. A modified quadratic cost problem and feedback stabilization of a linear system // IEEE Trans. Automat . Contr., AC - 22, 1977.

8Kwon W.N., Bruckstein A.N., Kailath T. Stabilizing state feedback design via the moving horizon method, Int. J. Contr., V. - 37. - 1983.

9Mayne D.Q., Michalska //. Receding horizon control of nonlinear systems // IEEE Trans. Automat. Contr, V.3, №7. - P. 814 - 824.

10R.Gabasov, F.M.KiriMova, E.A.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.

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

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

13Габасов Р. Применение позиционных решений кусочно-линейно-квадратичных задач для демпфирования и стабилизации маятника / Р. Габасов, Ф.М. Кириллова, А.В. Лубочкин // Известия РАН. Теория и системы управления. - 2002. - № 5. - С. 64 - 73.

14Габасов Р. Стабилизация в большом перевернутого маятника / Р. Габасов [и др.] //Известия РАН. Теория и системы управления. - 2003. - № 1. - С. 17 - 23.

15Габасов Р. Демпфирование и стабилизация маятника при больших начальныхвозмущениях / Р. Габасов, Ф.М. Кириллова, Е.А. Ружицкая // Известия РАН. Теория и системы управления. - 2001. - №1. - С. 29 - 38.

16Габасов Р., Кириллова Ф.М., Ружицкая Е.А. Стабилизация перевернутого маятника // Доклады НАН Беларуси. 2000. Т.44, т. С. 9 - 12.

17Furuta К., Yamakita М., Kobayashi S. Swing-up control of inverted pendulum using pseudo-state feedback // ,5. Systems and Control Eng. 1992. V.206, P. 263 - 269.

18Astrom K.J., Furuta K. Swinging up a. pendulum by energy control // Proc. 13th Triennial World Congress. San Francisco, USA. - P. 37 - 42.

19Конструктивные методы оптимизации / P. Габасов, Ф.М. Кириллова, О.П. Костюкова, А.В. Покатаем // Нелинейные задачи. - Мн.: Университетское, 1998. - Ч. 5. - 390 с.

20Gabasov R., Kirillova F.M. Real-time construction of optimal closable feedbacks 13th World Congress of IFAC International Federation of Automatic Control, V.D., San Francisca, CA, USA, 1996

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

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

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

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

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

26Лубочкин А.В. Реализация позиционного решения задачи оптимизации линейной системы по негладкому критерию качества / А.В. Лубочкин // Известие РАН. Теория и системы управления. - 1998. - № 1. - С. 64 - 70.

Приложение А


Код программы

# include <stdio.h>

# include <math.h>

# include <string.h>

# include <process.h>

# include <alloc.h>

# include <conio.h>

#define M 5

#define N 30

#define N1 N+1

#define TBEG 0.0

#define TETA 3.0

#define NX 2

#define pi 3.1415926536x[NX] = { pi, 0.0 },[N1],b[M],A[M][N1],[N1],dw[N1],u[N1],[N1], ur,J,ga;nop, Jop[M],StrZ,,NachRe=0,ProgrRe=0,,nnz;puN1NRo ( int*, int*, int*, double );puRo(int,int,int,double);formParam ( int, int, int, double );f_x (double, double*, // f1(t)=x2(t)* ); //f2(t)=sin(x1(t))- u(t)*cos(x1(t))

//или

//f1(t)=x2(t)

//f2(t)=a(x1(t))-u(t)*b(x1(t))FEHL(int, double,double*, void (f_p)(double, double*, double*),*, double*,*, double*,*, double*,*, double*);RKF45 (int, double, // Метод Рунге-Кутта-Фельдберга, double*,(f_p)(double, double*, double*), double*,*, double*,*, double*,*, double*,*, int*,*, int*,*, int*,*, double*,*, double*);prifl0 (int*, double*, double*);*d;*r;*ru; // Файл.dat : t, u(t)*rJ; // Файл.dat : t, J(t)*rx; // Файл.dat : t, x(t)main (void)

{char f_Name_u[13],// Имя файла.dat : t, u(t)_Name_J[13],// Имя файла.dat : t, J(t)_Name_x[13];// Имя файла.dat : t, x(t)tau,h,t,tout,ypx [NX],=0,f1x [NX],f2x [NX],x [NX],f4x [NX],x [NX],rerr=1.0e-12,= 0.0,savrex = 0,= 0;j,n1,n2,n3,Ostanov,inte,, jflx, kflx, nfex,, initx;

//Открытие файлов:( "\nFileName: " );( "%s", f_Name_u );( f_Name_x, f_Name_u );( f_Name_J, f_Name_u );( f_Name_u, "_u.dat" );( f_Name_J, "_J.dat" );( f_Name_x, "_x.dat" );((ru=fopen(f_Name_u,"wt")) == NULL )

{ perror ( f_Name_u );(-1);}((rJ = fopen ( f_Name_J, "wt")) == NULL )

{ perror ( f_Name_J );(-1);}((rx = fopen(f_Name_x, "wt")) == NULL )

{ perror ( f_Name_x );(-1);}

tau = TBEG;

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

fprintf ( rx, "%lf ", tau );( "\ntau = %lf",tau );( "x = ( " );( j=0; j<NX; j++ )

{fprintf ( rx, "%lf ",x[j] );( "%lf ", x[j] );}( rx, "\n" );( ")\n" );( "| x(1) - 3pi/4 | = %lf\n",fabs(x[0]-3.0*pi/4.0));( "| x(1) - 5pi/4 | = %lf\n\n",fabs(x[0]-5.0*pi/4.0) );

//Начальн. формирование некот. параметров конечномерн.задачи:= TETA/N; // шаг(длина перем.)

printf ( "h = %lf\n\n", h );= sqrt( 4.0/pi - 1.0 );= 1;( "n1 = " );( "%d", &n1 );( "n2 = " );( "%d", &n2 );( "n3 = " ); ( "%d", &n3 );

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

Ostanov = 0;

while ( !Ostanov )

{//Программн. реш. задачи с перем. структурой:

if ( ( StrZ==1 ) || ( StrZ==2 ) || ( StrZ==3 ) )(1)

{printf ( "\nСтруктура задачи:\n" );( "1 -- четыре области: [ 3pi/4, 5pi/4 ]\n" );

printf("[pi/2, 3pi/4]\n" );("[pi/4, pi/2]\n" );("[-pi/4,pi/4]\n\n"); ( "2 -- три области:

[pi/2, 3pi/4]\n" );

printf("[ pi/4, pi/2]\n" );("[-pi/4,pi/4]\n\n");

printf("3 -- две области:

[pi/4, pi/2]\n" );

printf("[-pi/4,pi/4]\n\n");

printf("4 -- одна область:

[-pi/4,pi/4]\n\n");

scanf ( "%d", &StrZ );( StrZ==1 )

{ m = M;;}( StrZ==2 )

{m = M-1;;}( StrZ==3 )

{ m = M-2;;}( StrZ==4 )

{ m = M-3;

break;}("\nНеправильно выбрана структура !!!\n\n" );}( StrZ==1 )

{printf("\nПараметры конечномерн.з-чи:\n" );

printf("Размеры: m=%d n1=%d n2=%d n3=%d N=%d\n", m, n1, n2, n3, N );

printf( "h = %lf\n", h );();();}( StrZ==2 )

{printf ( "\nПараметры конечномерн.з-чи:\n" );

printf ( "Размеры: m=%d n2=%d n3=%d N=%d\n",

m, n2, n3, N );( "h = %lf\n", h );();();}

if ( StrZ==3 )

{printf("\nПараметры конечномерн.з-чи:\n" );("Размеры: m=%d n3=%d N=%d\n", m, n3, N );( "h = %lf\n", h );();();}NRo ( &n1, &n2, &n3, h );

//Печ. в файл u(t), J(t), x(t) на [tau,tau+h]:( ru, "%lf %lf\n", tau, u[0] );(ru,"%lf %lf\n",tau+h, u[0]);(rJ,"%lf %lf\n",tau,J);= tau;( "tau = %lf\n", tau );( tout < tau+h-0.000001 )

{t = tout;( tau + h - tout > 0.1 )+= 0.1;= tau + h;( "\nt = %lf", t );("tout = %lf\n",tout );

//getchar();= 1;= 1;(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);( iflx==2 )= 0;

{printf ("x\n");(&iflx, &rerr, &aerr);}}( rx, "%lf ", tout );( "x = ( " );( j=0; j<NX; j++ )

{fprintf( rx, "%lf ",x[j]);( "%lf ", x[j] );}( rx, "\n" );( ") " );( StrZ==1 )

{printf ( "\n" );("|x(1)-3pi/4|=%lf\n",( x[0]-3.0*pi/4.0 ) );("|x(1)-5pi/4|=%lf\n",( x[0]-5.0*pi/4.0));}( StrZ==2 )("| x(1) - pi/2 |=%lf\n",fabs( x[0]-pi/2.0));( StrZ==3 )("| x(1) - pi/4 | = %lf\n",fabs( x[0]-pi/4.0));( StrZ==4 )("| x(1) | = %lf\n", fabs( x[0] )); }+= h;( "\n\n\n tau = %lf\n\n\n", tau );

while (1)

{printf ( "Остановить процесс ? ( 1 -- да, 0 -- нет )" );

scanf ( "%d", &Ostanov );((Ostanov==0)||(Ostanov==1));}

} // while( !Ostanov )();

} // mainpuN1NRo ( int *n1, int *n2, // Прогр.упр.задачи пер.стр.*n3, double h )

{double J0;i, nz1, nz2, nz3,,MinLeft,,Uvelich,;

nz1 = *n1;= *n2;= *n3;

//Программн. упр-ие мин. инт.:( ( StrZ==1 ) || ( StrZ==2 ) || ( StrZ==3 ) )

{//Нач. формир. парам. n1,n2,n3

if (( StrZ==1) && (nz1 < m-3))

{if ( !ProgrRe )

{printf ( "\nНачало прогр. решения:\n" );

ProgrRe = 1;= 0;= nz1;

} // for if ( !ProgrRe )( ProgrRe )

{if ( ii < nnz )

{ ii++;[0] = u[ii];( "\nu0[] =\n" );(i=0; i<N; i++)

{printf ( "%10.6lf", u[i] );} ( "\n\nПрограммное решение:\n" );

printf ( "\nu0[0] = u0[%d] = %10.6lf\n",ii,u[0]);();

} // for if ( ii < nnz1 )

if ( ii == nnz )

{ printf ( "\n\nКонец программного решения\n\n" );

getchar();= 0;}

} // for if ( ProgrRe )endfunc;

} // for if ( ( StrZ==1 ) && ( nz1 < m-3 ) )((StrZ==2) && (nz2 < m-2))

{if ( !ProgrRe )

{ printf ( "\nНачало программного решения:\n" );

ProgrRe = 1;= 0;= nz2;

} // for if ( !ProgrRe )( ProgrRe )

{ ii++;[0] = u[ii];( "\nu0[] =\n" );(i=0; i<N; i++)

{printf("%10.6lf", u[i]);}( "\n\nПрограммное решение:\n" );( "\nu0[0] = u0[%d] = %10.6lf\n",ii, u[0]);();

} // for if ( ii < nnz2 )

if ( ii == nnz )

{printf ( "\n\nКонец программного решения\n\n" );

getchar();= 0;}

} // for if ( ProgrRe )endfunc;

} // for if ((StrZ==2)&&(nz2< m-2))((StrZ==3) && (nz3 < m-1))

{if ( !ProgrRe )

{ printf ( "\nНачало программного решения:\n" );

ProgrRe = 1;= 0;= nz3;

} // for if ( !ProgrRe )( ProgrRe )

{if ( ii < nnz )

{ ii++;[0] = u[ii];( "\nu0[] =\n" );(i=0; i<N; i++)( "%10.6lf", u[i] );} ( "\n\nПрограммное решение:\n" );

printf ( "\nu0[0] = u0[%d] = %10.6lf\n",ii, u[0]);();

} // for if ( ii < nnz3 )

if ( ii == nnz )

{ printf ( "\n\nКонец программного решения\n\n" );

getchar();= 0;}

} // for if ( ProgrRe )endfunc;

} // for if ( ( StrZ==3 ) && ( nz3 < m-1 ) ):NachRe = 0;= 0;( nz1, nz2, nz3, h );( nop==m )

{for ( i=0; i<N; i++ )[i] = u[i];= 1.0/u[N];= 1;}

{if ( StrZ==1 )("При n1=%d n2=%d n3=%d решение не найдено!\n",nz1,nz2,nz3);( StrZ==2 )("При n2=%d n3=%d решение не найдено\n",nz2,nz3);( StrZ==3 )("При n3=%d решение не найдено !!!\n",nz3);}:if ( !NachRe )( StrZ==1 )

{ printf ( "n1 = " );( "%d", &nz1 );( "n2 = " );( "%d", &nz2 );( "n3 = " );( "%d", &nz3 );mm1;}( StrZ==2 )

{ printf ( "n2 = " );( "%d", &nz2 );( "n3 = " );( "%d", &nz3 );mm1;}( StrZ==3 )

{ printf ( "n3 = " );( "%d", &nz3 );

goto mm1;}


//Оптимизация парам. n1, n2,n3

if ( StrZ==1 )

{printf ( "n1 = %d n2 = %d n3 = %d\n", nz1, nz2, nz3 );("J = %.12lf\n\n\n",J0);

//getchar();= 0; // <--- ---> nz1(1) // <--- nz1

{if ( nz1 > 1 )

{while (1)

{printf("Продолжить уменьшение n1?( 1 -- да, 0 -- нет)");

scanf ( "%d", &Umensh );( ( Umensh==0 ) || ( Umensh==1 ) );}( Umensh==0 );-= 1;("Размеры: m=%d n1=%d n2=%d n3=%d N=%d\n",m,nz1,2,nz3,N);

getchar();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n1 = %d n2 = %d n3 = %d\n",nz1, nz2, nz3 );("J = %.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;= 1;}

{nz1 += 1;;}}// for if ( nop==m )

{nz1 += 1;;}}// for if ( nz1 > 1 );} // while(1)( !MinLeft ) // ---> nz1

{while(1)

{if ( nz1 < nz2-1 )

{while (1)

{ printf ( "Увеличить n1 ? ( 1 -- да, 0 -- нет )" );

scanf ( "%d", &Uvelich );

if ((Uvelich==0)||(Uvelich==1));}( Uvelich==0 );+= 1;("Размеры: m=%d n1=%d n2=%d n3=%d N=%d\n",m, nz1,, nz3, N);();

//Программн. упр-ие мин. инт.:= 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n1 = %d n2 = %d n3 = %d\n",nz1,nz2,nz3 );("J = %.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;}

{ nz1 -= 1;;}}// for if ( nop==m )

{ nz1 -= 1;;}}// for if ( nz1 < nz2-1 );

} // while(1)

} // if ( !MinLeft ) // end for <--- ---> nz1:MinLeft = 0; // <--- ---> nz2(1) // <--- nz2

{if ( nz2 > nz1+1 )

{while (1)

{printf("Продолжить уменьшение n2? (1 -- да, 0 -- нет )");

scanf ( "%d", &Umensh );

if ((Umensh==0)||(Umensh==1));}( Umensh==0 );-= 1;("Размеры: m=%d n1=%d n2=%d n3=%d N=%d\n", m, nz1, nz2, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n1 = %d n2 = %d n3 = %d\n",nz1, nz2, nz3 );("J = %.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;= 1;}

{ nz2 += 1;;}}//for if( nop==m )

{ nz2 += 1;;}}// for if (nz2>nz1+1);

} // while(1)( !MinLeft ) // ---> nz2

{while(1)

{if ( nz2 < nz3-1 )

{while (1)

{printf ( "Увеличить n2 ?

( 1 -- да, 0 -- нет )" );

scanf ( "%d", &Uvelich );

if ( ( Uvelich==0 ) || ( Uvelich==1 ) );}( Uvelich==0 );+= 1;("Размеры: m=%d n1=%d n2=%d n3=%d N=%d\n", m, nz1, nz2, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];("n1 = %d n2 = %d n3 = %d\n",nz1, nz2, nz3);("J = %.12lf\n\n\n",J); getchar();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;}

{ nz2 -= 1;;}}// for if ( nop==m )

{ nz2 -= 1;;}}// for if (nz2<nz3-1);

} // while(1)

} // if ( !MinLeft ) // end for <--- ---> nz2= 0; // <--- ---> nz3(1) // <--- nz3

{if ( nz3 > nz2+1 )

{while (1)

{printf("Продолжить уменьшение n3?(1--да,0--нет)");

if ((Umensh==0)||(Umensh==1));}( Umensh==0 );-= 1;("Размеры: m=%d n1=%d n2=%d n3=%d N=%d\n", m, nz1, nz2, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n1 = %d n2 = %d n3 = %d\n",nz1, nz2, nz3 );("J = %.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;= 1; }

{ nz3 += 1;;}}// for if ( nop==m )

{ nz3 += 1;;}}//for if (nz3>nz2+1);

} // while(1)( !MinLeft ) // ---> nz3

{while(1)

{if ( nz3 < N-1 )

{while (1)

{printf ( "Увеличить n3 ? ( 1 -- да, 0 -- нет )" );

scanf ( "%d", &Uvelich );

if ((Uvelich==0)||(Uvelich==1));}( Uvelich==0 );+= 1;("Размеры: m=%d n1=%d n2=%d n3=%d N=%d\n", m, nz1, nz2, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];("n1 = %d n2 = %d n3 = %d\n",nz1, nz2, nz3 );("J=%.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;}

{ nz3 -= 1;;}}// for if ( nop==m )

{ nz3 -= 1;;}}//for if (nz3< N-1);

} // while(1)

} // if ( !MinLeft ) // end for <--- ---> nz3( i=0; i<N; i++ )[i] = u0[i];= J0;( "n1 = %d n2 = %d n3 = %d\n", nz1, nz2, nz3 );("J = %.12lf\n\n\n",J);();= 0;("\nИзменить n1, n2, n3?( 1 -- да, 0 -- нет)");( "%d", &KorrektN );( KorrektN )

{ printf ( "n1 = " );( "%d", &nz1 );( "n2 = " );( "%d", &nz2 );( "n3 = " );( "%d", &nz3 );= 0;( nz1, nz2, nz3, h );( nop==m )

{for ( i=0; i<N; i++ )[i] = u[i];= 1.0/u[N];( "n1 = %d n2=%d n3 = %d\n",nz1, nz2, nz3 );("J= %.12lf\n\n\n",J0);();m1;}// for if ( nop==m )mm2;}

} // if ( StrZ==1 )( StrZ==2 )

{printf ( "n2 = %d n3 = %d\n", nz2, nz3 );("J=%.12lf\n\n\n",J0 );

//getchar();= 0; /* <--- ---> nz2 */ while(1) /* <--- nz2 */

{if ( nz2 > 1 )

{while (1)

{printf("Продолжить уменьшение n2 ? ( 1 -- да, 0 -- нет )");( "%d", &Umensh );

if ( ( Umensh==0 ) || ( Umensh==1 ) );}( Umensh==0 );-= 1;( "Размеры: m=%d n2=%d n3=%d N=%d\n",m, nz2, nz3, N );();

//Программн. упр-ие миним. инт.

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n2 = %d n3 = %d\n", nz2, nz3 );("J=%.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;= 1;}

{ nz2 += 1;;}}// for if ( nop==m )

{ nz2 += 1;;}}// for if (nz2 > 1);

} // while(1)( !MinLeft ) // ---> nz2

{while(1)

{if ( nz2 < nz3-1 )

{while (1)

{printf ( "Увеличить n2 ? ( 1 -- да, 0 -- нет )" );

scanf ( "%d", &Uvelich );

if ((Uvelich==0)||(Uvelich==1));}( Uvelich==0 );+= 1;( "Размеры: m=%d n2=%d n3=%d N=%d\n",m, nz2, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n2 = %d n3 = %d\n", nz2, nz3 );("J = %.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;}

{ nz2 -= 1;;}}// for if ( nop==m )

{ nz2 -= 1;;}}// for if (nz2<nz3-1);

} // while(1)

} // if ( !MinLeft ) end for <--- ---> nz2: MinLeft = 0; // <--- ---> nz3(1) // <--- nz3

{if ( nz3 > nz2+1 )

{while (1)

{printf("Продолжить уменьшение n3 ? ( 1 -- да,0 -- нет)");

scanf ( "%d", &Umensh );

if ((Umensh==0)||(Umensh==1));}( Umensh==0 );-= 1;( "Размеры: m=%d n2=%d n3=%d N=%d\n",m, nz2, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n2 = %d n3 = %d\n", nz2, nz3 );("J = %.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;= 1;

{nz3 += 1;;}}// for if ( nop==m )

{ nz3 += 1;;}}// for if (nz3>nz2+1);

} // while(1)( !MinLeft ) // ---> nz3

{while(1)

{if ( nz3 < N-1 )

{while (1)

{ printf ( "Увеличить n3 ? ( 1 -- да, 0 -- нет )" );

scanf ( "%d", &Uvelich );

if ((Uvelich==0)||(Uvelich==1));}( Uvelich==0 );+= 1;( "Размеры: m=%d n2=%d n3=%d N=%d\n", m, nz2, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n2 = %d n3 = %d\n", nz2, nz3 );("J= %.12lf\n\n\n",J );();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;}

{ nz3 -= 1;;}}// for if ( nop==m )

{ nz3 -= 1;;}}// for if ( nz3 < N-1 );

} // while(1)

} // if ( !MinLeft ) end for <--- ---> nz3( i=0; i<N; i++ )[i] = u0[i];= J0;( "n2 = %d n3 = %d\n", nz2, nz3 );("J=%.12lf\n\n\n",J);();= 0;( "\nИзменить n2, n3? ( 1 -- да, 0 -- нет )" );

scanf ( "%d", &KorrektN );

if ( KorrektN )

{printf ( "n2 = " );( "%d", &nz2 );( "n3 = " );( "%d", &nz3 );= 0;( nz1, nz2, nz3, h );( nop==m )

{for ( i=0; i<N; i++ )[i] = u[i];= 1.0/u[N];}( "n2 = %d n3 = %d\n", nz2, nz3 );("J = %.12lf\n\n\n",J0);();m2;}

} // if ( StrZ==2 )( StrZ==3 )

{printf ( "n3 = %d\n", nz3 );("J=%.12lf\n\n\n",J0 );

//getchar();: MinLeft = 0; // <--- ---> nz3(1) // <--- nz3

{if ( nz3 > 1 )

{while (1)

{printf("Продолжить уменьшение n3 ? ( 1 -- да, 0 -- нет )");( "%d", &Umensh );

if ((Umensh==0)||(Umensh==1));}( Umensh==0 );-= 1;( "Размеры: m=%d n3=%d N=%d\n",m, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n3 = %d\n", nz3 );("J=%.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;= 1;}

{nz3 += 1;;}}// for if ( nop==m )

{ nz3 += 1;;}}// for if ( nz3 > 1 );

} // while(1)( !MinLeft ) // ---> nz3

{while(1)

{if ( nz3 < N-1 )(1)

{ printf ( "Увеличить n3 ? ( 1 -- да, 0 -- нет )" );

scanf ( "%d", &Uvelich );

if ((Uvelich==0)||(Uvelich==1));}( Uvelich==0 );+= 1;("Размеры: m=%d n3=%d N=%d\n",m, nz3, N );();

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

nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];( "n3 = %d\n", nz3 );("J=%.12lf\n\n\n",J);();( J<J0 )

{for ( i=0; i<N; i++ )[i] = u[i];= J;}

{ nz3 -= 1;;}}// for if ( nop==m )

{ nz3 -= 1;;}}// for if ( nz3 < N-1 );

} // while(1)

} // if ( !MinLeft ) end for <--- ---> nz3( i=0; i<N; i++ )[i] = u0[i];= J0;( "n3 = %d\n", nz3 );("J = %.12lf\n\n\n",J);();= 0;( "\nИзменить n3? ( 1 -- да, 0 -- нет ) " );

scanf ( "%d", &KorrektN );

if ( KorrektN )

{ printf ( "n3 = " );( "%d", &nz3 );= 0;( nz1, nz2, nz3, h );( nop==m )

{for ( i=0; i<N; i++ )[i] = u[i];= 1.0/u[N];}( "n3 = %d\n", nz3 );("J=%.12lf\n\n\n",J0);();m3;}

} // if ( StrZ==3 )}// for if ((StrZ==1)||

(StrZ==2 ) || ( StrZ==3 ) )

{ nop = 0;( nz1, nz2, nz3, h );( nop==m )

{J = 1.0/u[N];("J=%.12lf\n",J );();}

{printf ("Нет решения !!!\n" );();}

} //end else for if (( StrZ==1 ) || ( StrZ==2 ) || ( StrZ==3 )):

*n1 = nz1-1;

*n2 = nz2-1;

*n3 = nz3-1;

} // puN1NRopuRo ( int n1, int n2, // Прогр.упр. мин.инт.n3, double h )

{ int i, j, k;

//Формирование параметров конечномерн. задачи:( n1, n2, n3, h );

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

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

if ( ( d = fopen ( "d.d", "wt" ) ) == NULL )

{ perror ( "d.d" );(-1);}

//Печать параметров задачи в файл(d,"%d %d\n", m, N1 ); // размеры m и N1

for ( i=0; i<m; i++ ) // матр.A, вект.b

{for ( j=0; j<N1; j++ )( d, "%12.6lf ", A[i][j] );( d, "\n" );( d, "%11.6lf %11.6lf\n", b[i], b[i] );}( j=0; j<N1; j++ ) // вект. dn, dw, c( d, "%4.1lf %4.1lf %11.6lf\n", dn[j], dw[j], c[j] );(d);( "lplbig.exe" );

//printf ( "system()\n");

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

if ( ( r = fopen ( "sol.d", "r" ) ) == NULL )

{perror ( "sol.d" );

exit (-1);}

fscanf ( r, "%d", &nop ); // Размер опоры

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

{ fscanf ( r, "%d", &k );[i] = k-1;}(j=0; j<N1; j++)

{ fscanf ( r, "%lf", &ur );[j] = ur;}(j=0; j<N; j++)

{ u[j] = u[j]/u[N];}( "\nu0[0] = %10.6lf\n", u[0] );(r);( "\nu0[] =\n" );( i=0; i<10; i++ )( "%7.1lf", u[i] );( "\n" );( i=10; i<20; i++ )( "%7.1lf", u[i] );( "\n" );( i=20; i<N; i++ )( "%7.1lf", u[i] );( "\n" );

//getchar();

} // puRoformParam ( int n1, int n2, // Форм.парам.конечном.з-чиn3, double h )

{ double t1_tn,t1_tw,t2_tn,_tw,t3_tn,t3_tw,_tn,t_tw,t2_t1,_t2,t_t3,,t1,t2,t3;i, j;( j=0; j<N; j++ )[j] = 0.0;[N] = 1.0;( StrZ==1 )

{t2_t1 = n2*h - n1*h;_t2 = n3*h - n2*h;_t3 = N*h - n3*h;= n1*h;

for ( j=0; j<n1; j++ ) // матрица A

{t1_tn = n1*h - j*h;_tw = n1*h - (j+1)*h;[0][j]=cos(t1_tw)-cos(t1_tn);[1][j]=sin(ga*t2_t1)*(sin( t1_tn)-sin(t1_tw))/ga;[2][j]=sinh(ga*t3_t2)*cos(ga*t2_t1)*(sin(t1_tn)-sin(t1_tw))/ga;[3][j]=sinh(t_t3)*cosh(ga*t3_t2)*cos(ga*t2_t1)*(sin(t1_tn)-sin(t1_tw));[4][j]=cosh(t_t3)*cosh(ga*t3_t2)*cos(ga*t2_t1)*(sin(t1_tn)- sin(t1_tw));}( j=n1; j<n2; j++ )

{t2_tn = n2*h - j*h;_tw = n2*h - (j+1)*h;[0][j]=0.0;[1][j]=cos(ga*t2_tw)-cos(ga*t2_tn);[2][j]=sinh(ga*t3_t2)*sin(ga*t2_tn)-sin(ga*t2_tw));[3][j]=ga*sinh(t_t3)*cosh(ga*t3_t2)*(sin(ga*t2_tn)-sin(ga*t2_tw));[4][j]=ga*cosh(t_t3)*cosh(ga*t3_t2)*(sin(ga*t2_tn)-(ga*t2_tw ));}( j=n2; j<n3; j++ )

{t3_tn = n3*h - j*h;_tw = n3*h - (j+1)*h;[0][j]=0.0;[1][j]=0.0;[2][j]=cosh(ga*t3_tw)-cosh(ga*t3_tn);[3][j]=ga *sinh(t_t3)*(sinh(ga*t3_tw-sinh(ga*t3_tn));[4][j]=ga*cosh(t_t3)*(sinh(ga*t3_tw)-sinh(ga*t3_tn));}( j=n3; j<N; j++ )

{t_tn = N*h - j*h;_tw = N*h - (j+1)*h;[0][j]=0.0;[1][j]=0.0;[2][j]=0.0;[3][j]=cosh(t_tw)-cosh(t_tn);[4][j]=sinh(t_tw -sinh(t_tn);}[0][N]=-(3.0*pi/4.0-cos(t1)* x[0]-sin(t1)*x[1]-pi*(1.0-cos(t1)));[1][N]=-(pi/2.0-3.0*pi*cos(ga*t2_t1)/4.0+sin(ga*t2_t1)*sin(t1)*[0]/ga-sin(ga*t2_t1)*cos(t1)*x[1]/ga-pi*sin(ga*t2_t1)* sin(t1)/ga-(3.0-pi/2.0)*(1.0-cos(ga*t2_t1))/(ga*ga));[2][N]=-(pi/4.0-pi*cosh(ga*t3_t2)/2.0+3.0*pi*sinh(ga*t3_t2 )* sin(ga*t2_t1)/4.0+sinh(ga*t3_t2)*cos(ga*t2_t1)*sin(t1)*x[0]/ga-sinh(ga*t3_t2)*cos(ga*t2_t1)*(t1)*x[1]/ga-pi*(ga*t3_t2)* cos(ga*t2_t1)*sin(t1)/ga-(3.0-pi/2.0)*sinh(ga*t3_t2)*(ga*t2_t1 )/(ga*ga));[3][N]=-(-pi*cosh(t_t3)/4.0-pi*ga*sinh(t_t3)*sinh(ga*t3_t2)/2.0+3.0*pi*ga*sinh(t_t3)*(ga*t3_t2)*sin(ga*t2_t1)/4.0+sinh(t_t3)*cosh(ga*t3_t2)*cos(ga*t2_t1)*sin(t1)*x[0]-sinh(t_t3)*cosh(ga*t3_t2)*(ga*t2_t1)*cos(t1)*x[1]-pi*sinh(t_t3)*cosh(ga*t3_t2)*(ga*t2_t1)*sin(t1)-(3.0-pi/2.0)*sinh(t_t3)*cosh(ga*_t2)*sin(ga*t2_t1)/ga-(pi/2.0-1.0)*sinh(t_t3)*(ga*t3_t2)/ga );[4][N]=-(-pi*sinh(t_t3)/4.0-pi*ga*cosh(t_t3)*sinh(ga*t3_t2)/2.0+3.0*pi*ga*cosh(t_t3)*(ga*t3_t2)*sin(ga*t2_t1)/4.0+cosh(t_t3)*cosh(ga * t3_t2)*(ga*t2_t1)*sin(t1)*x[0]-cosh(t_t3)*cosh(ga*t3_t2)*(ga*t2_t1)*cos(t1)*x[1]-pi*cosh(t_t3)*cosh(ga*t3_t2)*os(ga*t2_t1)*sin(t1)-(3.0-pi/2.0)*cosh(t_t3)*(ga*t3_t2)*sin(ga*t2_t1)/ga-(pi/2.0-1.0)*cosh(t_t3)*( ga*t3_t2)/ga);}( StrZ==2 )

{t3_t2 = n3*h - n2*h;_t3 = N*h - n3*h;= n2*h;( j=0; j<n2; j++ )

{ t2_tn = n2*h - j*h;_tw = n2*h - (j+1)*h;[0][j]=os(ga*t2_tw)-cos(ga*t2_tn);[1][j]=sinh(ga*t3_t2)*(sin(ga*t2_tn)-sin(ga*t2_tw));[2][j]=ga*sinh(t_t3)*cosh(ga *t3_t2)*sin(ga*t2_tn)- sin(ga*t2_tw));[3][j]=ga*cosh(t_t3)*cosh(ga*3_t2)*(sin(ga*t2_tn)-(ga*t2_tw)); }( j=n2; j<n3; j++ )

{t3_tn = n3*h - j*h;_tw = n3*h - (j+1)*h;[0][j]=0.0;[1][j]=cosh(ga*t3_tw)-cosh(ga* t3_tn);[2][j]=ga*sinh(t_t3)*(sinh(ga*t3_tw)-sinh(ga*t3_tn));[3][j]=ga*cosh(t_t3)*(sinh(ga*t3_tw)-sinh(ga*t3_tn));}( j=n3; j<N; j++ )

{t_tn = N*h - j*h;_tw = N*h - (j+1)*h;[0][j]=0.0;[1][j]=0.0;[2][j]=cosh(t_tw)-cosh(t_tn);[3][j]=sinh(t_tw)-sinh(t_tn);}

// вектор b[0][N]=-(pi/2.0-cos(ga*t2)*x[0]-sin(ga*t2)*x[1]/ga-(3.0-pi/2.0)*

(1.0-cos(ga*t2))/(ga*ga));[1][N]=-(pi/4.0-pi*cosh(ga*t3_t2)/2.0+sinh(ga*t3_t2)*sin(ga*t2)* x[0]-sinh(ga*t3_t2)*cos(ga*t2)*x[1]/ga-(3.0-pi/2.0*sinh(ga*t3_t2) *sin(ga*t2)/(ga*ga)-(pi/2.0-1.0)*(cosh(ga*t3_t2)-1.0)/(ga*ga));[2][N]=-(-pi*cosh(t_t3)/4.0-pi*ga*sinh(t_t3)*sinh(ga*t3_t2)/2.0+ga*sinh(t_t3)*(ga*t3_t2)*sin(ga*t2)*x[0]-sinh(t_t3)*cosh(ga* t3_t2)*cos(ga*t2)*x[1]-(3.0-pi/2.0)*sinh(t_t3)*cosh(ga*t3_t2)*(ga*t2)/ga-(pi/2.0-1.0)*sinh(t_t3)*sinh(ga*t3_t2)/ga);[3][N]=-(-pi*sinh( t_t3)/4.0-pi*ga*cosh(t_t3)*(ga*t3_t2)/2.0+ga*cosh(t_t3)*cosh(ga*t3_t2)*sin(ga*t2)*x[0]-cosh(t_t3)*cosh(ga*t3_t2)*(ga*t2)*x[1]-(3.0-pi/2.0)*cosh(t_t3)*(ga *t3_t2)*sin(ga*t2)/ga-(pi/2.0-1.0)*cosh(t_t3)*(ga*t3_t2)/ga;}( StrZ==3 )

{t_t3 = N*h - n3*h;= n3*h;( j=0; j<n3; j++ )

{t3_tn = n3*h - j*h;_tw = n3*h - (j+1)*h;[0][j]=cosh(ga*t3_tw)-cosh(ga*t3_tn);[1][j]=ga*sinh(t_t3)*sinh(ga*t3_tw)-sinh(ga*t3_tn));[2][j]=ga*cosh(t_t3)*(sinh(ga*t3_tw)-sinh(ga*t3_tn));}( j=n3; j<N; j++ )

{t_tn = N*h - j*h;_tw = N*h - (j+1)*h;[0][j]=0.0;[1][j]=cosh(t_tw)-cosh(t_tn);[2][j]=sinh(t_tw)-sinh(t_tn);}

// вектор b[0][N]=-(pi/4.0-cosh(ga*t3)*x[0]-sinh(ga*t3)*x[1]/ga- (pi/2.0-1.0)*(cosh(ga*t3)-1.0)/(ga*ga));[1][N]=-(-pi*cosh(t_t3)/4.0-ga*sinh(t_t3)*sinh(ga*t3)*x[0]-(t_t3)*cosh(ga*t3)*x[1]-(pi/2.0-1.0)*sinh(t_t3)*(ga*t3)/ ga);[2][N]=-(-pi*sinh(t_t3)/4.0-ga*cosh(t_t3)*sinh(ga*t3)*x[0]-(t_t3)*cosh(ga*t3)*x[1]-(pi/2.0-1.0)*cosh(t_t3)*(ga*t3)/ga ); }( StrZ==4 )

{te = N*h;( j=0; j<N; j++ )

{ t_tn = N*h - j*h;_tw = N*h - (j+1)*h;[0][j]=cosh(t_tw)-cosh(t_tn);[1][j]=sinh(t_tw)-sinh(t_tn);}[0][N]=-(-cosh(te)*x[0]-sinh(te)*x[1]);[1][N]=-(-sinh(te)*x[0]-cosh(te)*x[1]);}( i=0; i<m; i++ )[i] = 0.0;( j=0; j<N; j++ )

{ dn[j] = -1.0; dw[j] = 1.0;}[N] = -100000000.0;[N] = 100000000.0;

} // formParam()f_x ( double t, // Выч. f1(t) = x2(t)y[], // f2(t) = -sin(x1(t))+u(t)f[] )

{double ax, bx;

/*( ( y[0] > -pi/4.0 ) && ( y[0] <= pi/4.0 ) )

{ ax = y[0];bx = 1.0; }( ( y[0] > pi/4.0 ) && ( y[0] <= pi/2.0 ) )

{ ax = ga*ga * y[0] + pi/2.0 - 1.0; bx = ga*ga;}( ( y[0] > 3.0*pi/4.0 ) && ( y[0] <= 5.0*pi/4.0 ) )

{ ax=-y[0]+pi; bx=-1.0;}[0] = y[1];[1] = ax - bx * u[0];

*/[0] = y[1];[1] = sin( y[0] ) - cos( y[0] ) * u[0];

} // f_x()

# define True 1

# define False 0FEHL(int neqn, double t,y[], void (f_p)(double, double*, double*),*hc, double yp[],f1[], double f2[],f3[], double f4[],f5[], double s[])

{ double h;ch,rab;k;(*fun)(double, double*, double*);=*hc;= f_p;= h/4.0;(k=0; k<neqn; k++)[k] = y[k] + ch*yp[k];= t + ch;(rab,f5,f1);= 3.0 * h / 32.0;(k=0; k<neqn; k++)[k] = y[k] + ch*(yp[k] + 3.0 * f1[k]);= t + 3.0 * h / 8.0;(rab,f5,f2);= h / 2197.0;(k=0; k<neqn; k++)[k] = y[k] + ch * (1932.0 * yp[k] + (7296.0 * f2[k] - 7200.0 * f1[k]));= t + 12.0 * h / 13.0;(rab,f5,f3);= h / 4104.0;(k=0; k<neqn; k++)[k] = y[k] + ch * ((8341.0 * yp[k] - 845.0 * f3[k]) + (29440.0 * f2[k] - 32832.0 * f1[k]));= t + h;(rab,f5,f4);= h / 20520.0;(k=0; k<neqn; k++)[k] = y[k] + ch * ((-6080.0 * yp[k] + (9295.0 * f3[k] - 5643.0 * f4[k])) + (41040.0 * f1[k] - 28352.0 * f2[k]));= t + h / 2.0;(rab,f1,f5);

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

ch = h / 7618050.0;(k=0; k<neqn; k++)[k] = y[k] + ch *((902880.0 * yp[k] + (3855735.0 * f3[k] - 1371249.0 * f4[k]))+(3953664.0 * f2[k] + 277020.0 * f5[k]));

*hc=h;;}RKF45 (int neqn, double t,tout, double y[],(f_p)(double, double*, double*), double *hc,yp[], double f1[],f2[], double f3[],f4[], double f5[],*iflagc, int *nfec,*kopc, int *initc,*jflagc, int *kflagc,*relerrc, double *abserrc,double *savrec, double *savaec)

{ int iflag,nfe,kop,init,jflag,;h,savre,savae,, abserr;maxnfe = 30000;remin = 1.0e-15;hfaild, output;a,ae,dt,ee,eeoet,,et,hmin,eps,u26,rer,s,scale,tol,toln,ypk,rab1;k, mflag;(*fun)(double, double*, double*);= f_p;=*iflagc; nfe=*nfec; kop=*kopc; init=*initc; jflag=*jflagc;=*kflagc;=*hc; savre=*savrec; savae=*savaec;= *relerrc;= *abserrc;

// Проверить входные параметры(neqn < 1) goto m10;(relerr< 0.0||abserr<0.0)m10;=abs(iflag);(mflag == 0 || mflag > 8)m10;

if (mflag != 1) goto m20;

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

eps = 1.0;(eps+1.0 > 1.0)= eps/2.0;26 = 26.0*eps;

goto m50;

// Ошибки во входной информации: iflag = 8;mexit;

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

m20: if (t==tout && kflag != 3)m10;(mflag != 2) goto m25;(kflag == 3 || init==0)m45;(kflag == 4) goto m40;(kflag == 5 && abserr == 0.0) goto m30;(kflag == 6 && relerr <= savre && abserr <= savae)m30;m50;:if (iflag == 3)goto m45;(iflag == 4) goto m40;(iflag==5 && abserr > 0.0) goto m45;: goto mexit;: nfe = 0;(mflag == 2)m50;: iflag = jflag;(kflag == 3)= abs(iflag);: jflag = iflag;= 0;= relerr;= abserr;= 2.0*eps+remin;(relerr >= rer)m55;

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

relerr = rer;= 3;kflag = 3;mexit;: dt = tout-t;(mflag==1) goto m60;(init==0) goto m65;m80;: init = 0;= 0;= t;(a,y,yp);= 1;(t != tout)m65;= 2;mexit;: init = 1;= fabs(dt);= 0.0;(k=0; k<neqn; k++)

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

{ toln=tol;=fabs(yp[k]);=h*h*h*h*h;(ypk*rab1 > tol)=exp(0.2*log(tol/ypk));}}(toln <= 0.0) h=0.0;=u26*fabs(t);(fabs(t) < fabs(dt))=u26*fabs(dt);(h < rab1) h=rab1;=-2;(iflag > 0)=2;: if (dt < 0.0) h=-h;(fabs(h) >= 2.0*fabs(dt))=kop+1;(kop != 100 ) goto m85;

// Много выходов=0;=7;mexit;: if (fabs(dt) > u26*fabs(t))

goto m95;

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

for (k=0; k<neqn; k++)[k] = y[k] + dt * yp[k];=tout;(a,y,yp);=nfe+1;m300;: output=False;=2.0/relerr;=scale*abserr;

// Пошаговое интегрирование: hfaild=False;

hmin=u26*fabs(t);=tout-t;(fabs(dt) >= 2.0*fabs(h))m200;(fabs(dt) > fabs(h))m150;=True;=dt;

goto m200;: h=0.5*dt;

// Внутренний одношаговый интегратор

m200: if(nfe <= maxnfe)m220;

iflag=4;=4;mexit;

// Продолжить приближенное решение на один шаг длины h: FEHL(neqn, t, y, f_p, &h, yp, f1, f2, f3, f4, f5, f1);

nfe=nfe+5;=0.0;(k=0; k<neqn; k++)

{ et=fabs(y[k])+fabs(f1[k])+ae;

if (et > 0.0)m240;

// Неправильная граница погрешности

iflag=5;mexit;: ee=fabs((-2090.0*yp[k]+(21970.0*f3[k]-15048.0*f4[k]))+

(22528.0*f2[k]-27360.0*f5[k]));=ee/et;(eeoet < rab1)=rab1;}=fabs(h)*eeoet*scale/752400.0;

if (esttol <= 1.0)m260;

// Неудачный шаг; уменьшить его величину

hfaild=True;=False;=0.1;(esttol < 59049.0)=0.9/(exp(0.2*log(esttol)));=s*h;(fabs(h) > hmin)

goto m200;

// Неудача даже при минимальном шаге=6;=6;mexit;

// Успешный шаг

m260: t=t+h;(k=0; k<neqn; k++)[k]=f1[k];=t;(a,y,yp);=nfe+1;=5.0;(esttol > 1.889568e-4)=0.9/(exp(0.2*log(esttol)));(hfaild)

{ if (s > 1.0)=1.0; }=s*fabs(h);(rab1 < hmin)=hmin;(h > 0.0)=fabs(rab1);=-fabs(rab1);(output)m300;(iflag > 0)m100;

// Интегрирование успешно завершено=-2;mexit;: iflag=2;: *iflagc=iflag; *nfec=nfe; *kopc=kop; *initc=init; *jflagc=jflag;

*kflagc=kflag;

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

*relerrc = relerr;

*abserrc = abserr;;}prifl0 ( int *iflc,*rerrc,*aerrc )

{int ifl;rerr, aerr;= *iflc;= *rerrc;= *aerrc;(ifl)

{case 3: /*printf ("ifl=%d\n rerr=%.12lf aerr=%.12lf\n",ifl,,aerr);*/;4: printf ("много fp()\n ifl=%d\n", ifl);;5: aerr = 1.0e-9;

/*printf("ifl=%d\n rerr=%.12lf aerr=%.12lf\n",ifl,rerr,aerr);

*/;6: rerr *= 10.0;("ifl=%d\n rerr=%.12lf aerr=%.12lf\n", ifl, rerr, aerr)= 2;;7: printf ("много вых. точек\n ifl=%d\n", ifl);= 2;;8: printf ("ifl=%d\n", ifl);; }

*iflc = ifl;

*rerrc = rerr;

*aerrc = aerr;

//getchar();

}


МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ Учреждение образования «Гомельский государственный университет имени Франциска Скорины» Матем

Больше работ по теме:

КОНТАКТНЫЙ EMAIL: [email protected]

Скачать реферат © 2017 | Пользовательское соглашение

Скачать      Реферат

ПРОФЕССИОНАЛЬНАЯ ПОМОЩЬ СТУДЕНТАМ