close

Вход

Забыли?

вход по аккаунту

?

Narbut Labor rabot

код для вставкиСкачать
Л. К. НАРБУТ, О. В. БУКУНОВА
ЛАБОРАТОРНЫЕ РАБОТЫ ПО ТЕОРИИ
УПРАВЛЕНИЯ
Министерство образования и науки
Российской Федерации
Санкт-Петербургский государственный
архитектурно-строительный университет
Л. К. НАРБУТ, О. В. БУКУНОВА
ЛАБОРАТОРНЫЕ РАБОТЫ ПО ТЕОРИИ
УПРАВЛЕНИЯ
Учебное пособие
Санкт-Петербург
2013
1
УДК 517.926; 517.977
Предисловие
Рецензенты: д-р физ.-мат. наук, профессор Б. Г. Вагер (СПбГАСУ);
д-р физ.-мат. наук Е. Л. Генихович (Главная геофизическая обсерватория им. Воейкова)
Нарбут, Л. К.
Лабораторные работы по теории управления: учеб. пособие /
Л. К. Нарбут, О. В. Букунова; СПбГАСУ. – СПб., 2013. – 72 с.
ISBN 978-5-9227-0453-3
Пособие содержит четыре лабораторные работы. Первые две посвящены решению задачи Коши для системы двух линейных дифференциальных
уравнений 1-го порядка с постоянными коэффициентами. В третьей лабораторной работе производится построение траектории для линейной системы
двух уравнений, в четвертой – приближенное построение программного
управления, оптимального в смысле быстродействия, для системы двух линейных дифференциальных уравнений.
Описание лабораторной работы проводится по следующему плану:
 необходимые теоретические сведения;
 выполнение лабораторной работы;
 варианты для выполнения лабораторной работы.
Учебное пособие рассчитано на студентов, обучающихся по специальности «Прикладная математика».
Предмет «Теория управления» изучается студентами четвертого курса специальности «Прикладная математика» СПбГАСУ. Курс
состоит из двух частей. В первой части курса рассматриваются задачи, связанные с теорией устойчивости движения, с управляемостью и наблюдаемостью объектов, движение которых описывается
системами дифференциальных уравнений. Вторая часть курса посвящена решению задач, связанных с принципом максимума
Л. С. Понтрягина. В предлагаемое учебное пособие вошли лабораторные работы, выполняемые студентами при изучении второй части курса. Все работы выполняются на компьютере.
Авторы выражают глубокую благодарность профессору, заведующему кафедрой «Исследования операций» СПбГУ И. В. Романовскому за полезные консультации по программированию на языке Паскаль с использованием графического режима.
Авторы признательны рецензентам пособия д-ру физ.-мат. наук, профессору СПбГАСУ Б. Г. Вагеру и д-ру физ.-мат. наук, с. н. с.
Главной геофизической обсерватории им. Воейкова Е. Л. Гениховичу за ценные советы и полезные замечания.
Табл. 3. Ил. 10. Библиогр.: 7 назв.
Рекомендовано Редакционно-издательским советом СПбГАСУ в качестве учебного пособия.
ISBN 978-5-9227-0453-3
 Л. К. Нарбут, О. В. Букунова, 2013
 Санкт-Петербургский государственный
архитектурно-строительный университет, 2013
2
3
Введение
ЛАБОРАТОРНАЯ РАБОТА № 1
Предлагаемое пособие содержит четыре лабораторные работы.
Первая и вторая лабораторные работы посвящены решению задачи
Коши для системы двух линейных дифференциальных уравнений
1-го порядка с постоянными коэффициентами. В первой работе рассмотрена система однородных уравнений, а во второй – простейшая
система неоднородных уравнений. Обе работы выполняются сначала в системе MathCad, а затем на языке Паскаль.
Для выполнения третьей и четвертой работ используется язык
программирования Паскаль. В третьей лабораторной работе производится построение траектории для линейной системы двух уравнений. В четвертой лабораторной работе производится приближенное
построение программного управления, оптимального в смысле быстродействия, для системы двух линейных дифференциальных
уравнений.
Для всех лабораторных работ имеются индивидуальные задания (варианты 1–20). В конце книги есть список рекомендуемой литературы.
Решение задачи Коши для системы двух линейных
однородных дифференциальных уравнений 1-го порядка
с постоянными коэффициентами
4
Необходимые теоретические сведения
Как известно [1–3], уравнение движения управляемого объекта
в векторной форме имеет вид
(1.1)
x  f ( x, u),
где x – n -мерный вектор состояния объекта; u – r -мерный вектор
управления.
Если система линейна, то уравнение (1.1) принимает вид
(1.2)
x  Ax  Bu,
где А – матрица размера ( n  n ), В – матрица размера ( n  r ). В случае скалярного управления ( r  1) система (1.1) представляется в
виде
x  Ax  bu,
(1.3)
где b – n -мерный вектор.
В лабораторных работах будут рассматриваться линейные системы двух дифференциальных уравнений ( n  2 ) с одним управляющим параметром ( r  1). Они будут иметь вид (1.3) или
в скалярной форме
 x1  a11 x1  a12 x2  b1u,
(1.4)

 x 2  a21 x1  a22 x2  b2u.
Здесь aij – элементы матрицы А, bi – координаты вектора b
( i, j  1,2 ).
В лабораторных работах № 1 и 2 будут рассмотрены более
простые, по сравнению с (1.3), системы дифференциальных уравнений, не содержащие управление u . В лабораторной работе № 1 рассматривается система, векторная запись которой имеет вид
(1.5)
x  Ax
(линейная однородная система).
В лабораторной работе № 2 будет рассмотрена простейшая
линейная неоднородная система, также не содержащая управляющий параметр u , ее векторная запись
(1.6)
x  Ax  v,
где v – постоянный вектор размерности n  2 .
5
Начиная с лабораторной работы № 3, будут рассматриваться
системы, содержащие управление u .
Итак, в лабораторной работе № 1 рассматривается система
дифференциальных уравнений (1.4):
 x1  a11 x1  a12 x2 ,
(1.7)

 x 2  a21 x1  a22 x2 .
Требуется найти частное решение этой системы, удовлетворяющее начальным условиям:
(1.8)
x1 (0)  x10 , x2 (0)  x20 ,
т. е. найти решение задачи Коши.
Для решения этой задачи можно использовать известные методы: метод исключения, метод Эйлера, операционный метод, рассмотренные в [5, 6]. Эти методы удобны для решения конкретных
числовых примеров. При работе на компьютере удобнее использовать точные формулы, рассмотренные в [4, 6]. Остановимся подробнее на этих формулах. Задачи (1.7), (1.8) в векторной форме
имеет вид:
найти частное решение векторного уравнения
(1.9)
x  Ax ,
удовлетворяющее начальному условию
x0  x 0 .
(1.10)
Как известно, решение рассматриваемой задачи Коши представляется в виде
xt   e At x 0 .
(1.11)
At
Матричная экспонента e может быть найдена по одной из
трех формул, рассмотренных далее.
Случай 1. Собственные значения матрицы А вещественны
и различны, λ1  λ 2 :
A  2E
A  1E
.
(1.12)
e At  e 1t
 e  2t
1   2
 2  1
Случай 2. Собственные значения матрицы А совпадают,
λ1  λ 2 :
(1.13)
e At  e λ 0t E  A  λ 0 E t .
Случай 3. Собственные значения матрицы А комплексносопряженные, λ1, 2  α  βi :
xx2 t   2e  3t .
Обычно эти формулы студенты получают при изучении первой части курса «Теория управления», решая эту же задачу аналитически, например операционным методом. Эти формулы используются для сравнения результатов аналитического и компьютерного
решений. Такое сравнение позволяет проверить правильность работы программы.
6
7
A  E


(1.14)
e At  e  t  E cos  t 
sin t .



Замечание 1.1. Вывод формул, аналогичных формулам (1.12)–
(1.14), проводится в [4, 6].
Выполнение лабораторной работы № 1
Вариант 0
1. Постановка задачи
Дана система дифференциальных уравнений
 x1  2 x1  x2 ,

 x2  3x2 .
Требуется найти частное решение системы, удовлетворяющее
начальным условиям
x1 0   2, x 2 0   2 .
2. Решение в системе MathCad
Для решения задачи нам потребуются следующие данные:
1
 2
а) A  
 – матрица коэффициентов правой части сис 0  3
темы;
 2
б) x0    – вектор начальных условий;
 2
в) T – верхняя граница промежутка интегрирования, например, при T  1 аргумент t  0;1 ;
г) h – шаг аргумента t , например, при h  0.1 аргумент t принимает значения t  0 ; t  0.1 ; t  0.2 и так далее;
д) формулы составляющих искомого частного решения;
xx1 t   4e 2t  2e 3t ,
Текст программы
ORIGIN : 1
Единичная матрица
 1 0
E : 

 0 1
Матрица коэффициентов правой части системы
 2 1 
A : 

 0  3
Начальные условия
 2
x0 :  
 2
Верхняя граница промежутка интегрирования T  1
Шаг аргумента h : 0.1
T
Число шагов N : round 
N  10
h
Собственные значения матрицы A
  2
λ : eigenvalsA 
λ   
  3
 : Re1 
 : Im1 
  2
0
Нахождение матричной экспоненты:
λ2t
e 1 t
EXP1t  :
A  λ 2 E   e
A  λ1E 
1   2
 2  1
EXP 2t  : e λ1 t E  A  λ1E  t 


A  E
EXP 3t  : e  t  Ecos  t  
sin t 
β


EXPAt  : if λ1  λ 2 , EXP 2t , if β  0, EXP1t , EXP 3t 
Начальные значения искомых функций:
X11 : x 01
X 21 : x 0 2
i : 1.. N  1
Значения аргумента t i : h i  1
Искомое частное решение:
 X1i 

 : EXPAt i  x 0
 X 2i 
8
Решение, полученное аналитически:
XX1i : 4e  2t i  2e  3t i
XX2 i : 2e  3ti .
ti =
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
X1i =
2
1.793
1.584
1.382
1.195
1.025
0.874
0.741
0.626
0.527
0.442
X2i =
2
1.482
1.098
0.813
0.602
0.446
0.331
0.245
0.181
0.134
0.1
XX1i =
2
1.793
1.584
1.382
1.195
1.025
0.874
0.741
0.626
0.527
0.442
XX2i =
2
1.482
1.098
0.813
0.602
0.446
0.331
0.245
0.181
0.134
0.1
3. Решение на языке Паскаль
Для выполнения лабораторной работы № 1 на языке Паскаль
используется приводимая далее программа lr1, предназначенная для
решения задачи Коши для системы двух линейных однородных
уравнений с постоянными коэффициентами.
Решение задачи осуществляется с использованием точных
формул. Так как в языке Паскаль нет стандартной процедуры нахождения собственных значений матрицы, аналогичной функции
eigenvals в системе MathCad, а также нет операций с комплексными
числами и операций с матричными переменными, то все используемые при решении задачи формулы программируются с помощью
обычных средств языка Паскаль.
Для выполнения операций с матрицами (умножение матрицы
на вектор, сложение матриц, вычисление определителя и т. п.) составлены соответствующие функции и процедуры. Так же, как
в рассмотренной ранее программе, составленной в системе
MathCad, решение задачи находится с использованием формулы
(1.11), а также формул (1.12)–(1.14), по которым матричная экспонента находится в случаях вещественных различных собственных
значений матрицы А, вещественных совпадающих и комплексносопряженных собственных значений матрицы А.
9
Замечание 1.2. Конечно, приведенная ранее программа в системе MathCad значительно проще следующей далее программы на
языке Паскаль. Однако при решении задач оптимального управления, которые будут рассмотрены в лабораторных работах № 3, 4,
удобнее использовать язык программирования Паскаль, позволяющий программировать задачи со сложными алгоритмами и обладающий хорошими графическими возможностями. В основу этих программ будут положены приводимые далее программы, составленные
на языке Паскаль, для выполнения лабораторных работ № 1 и 2.
Исходные данные
массива xx с составляющими xx1, xx2 соответственно (см. далее
текст программы, строки 96–101).
Замечание 1.3. Исходные данные x min, x max, y min, y max
могут быть уточнены после вывода таблицы результатов (см. далее
результаты решения). В приведенных вариантах 1–20 для выполнения лабораторных работ обычно бывает достаточным задать значения x min  y min  3, x max  y max  3 .
Результаты решения
I. Результаты, связанные с вычислениями
Исходные данные задаются в разделе описаний констант
(см. далее текст программы, строки 4–8):
t 0 – нижняя граница промежутка интегрирования;
tt – верхняя граница промежутка интегрирования;
dt1 – шаг аргумента t при выводе результатов решения;
dt 2 – шаг аргумента t при построении графика;
x min, x max – нижняя и верхняя границы изменения переменной x1 (если нижняя граница неотрицательна, то x min задается равным некоторому отрицательному числу, например, x min  0.5 ; если
верхняя граница неположительна, то x max задается равным некоторому положительному числу, например, x max  0.5 );
y min , y max – нижняя и верхняя границы изменения переменной x 2 (в случае постоянства знака переменной x2 задается, например, y min  0.5 или y max  0.5 );
a – двумерный массив элементов матрицы A (матрицы коэффициентов при неизвестных в правой части системы);
x 0 – одномерный массив, вектор-столбец начальных условий;
grPath – константа строкового типа, путь к библиотеке графических подпрограмм (третий параметр в обращении к стандартной
процедуре InitGraph, в приведенном тексте grPath = ‘d : \ turbo’);
процедура CoshyResult введена для проверки правильности работы программы, задает точные формулы составляющих x1 (t ) ,
x2 (t ) частного решения, полученные заранее, без использования
данной программы, например, операционным методом; формулы
для составляющих x1 (t ) , x2 (t ) задаются с помощью одномерного
1. Сначала выводятся исходные данные: матрица А (массив a)
и начальный вектор x 0 (массив x0 ), см. строки 152–158.
2. Затем выводится число rr и собственные значения матрицы
А, см. строки 160–165:
а) число rr :
если собственные значения 1 и  2 матрицы А вещественны
и различны, то rr  1;
если собственные значения совпадают, то rr  2 ;
если собственные значения матрицы А комплексносопряженные, то rr  3 ;
б) собственные значения матрицы А:
если собственные значения 1 и  2 матрицы А вещественны,
то выводится результат
ll 1  1 ,
ll 2   2 ;
если собственные значения матрицы А комплексносопряженные, 1, 2     i .
3. После этого выводится таблица, состоящая из пяти столбцов, см. строки 167–168, 174:
 в 1-м столбце выводятся значения переменной t от t 0 до tt
с шагом dt1;
 во 2-м и 3-м столбцах выводятся соответствующие значения
составляющих x1 и x2 массива x – значения составляющих x1 t 
и x2 t  частного решения (они получаются в результате выполнения
процедуры Coshy);
10
11
 в 4-м и 5-м столбцах выводятся для контроля значения составляющих xx1 и xx2 массива xx – значения составляющих частного решения, формулы для которых задаются в процедуре
CoshyResult.
II. Результаты, связанные с построением графика
Строятся оси координат, на которых указываются значения
x min, x max, y min, y max, см. строки 182, 132–148;
строится траектория движения объекта из точки с координатами x10 , x20 по закону x1 t , x2 t  путем построения отрезков прямых, соединяющих точки с координатами x1, x2, построенные
последовательно в моменты времени из промежутка t 0; tt  с шагом
dt 2 , см. строки 185–191.
В программе имеются комментарии, которые записаны в приведенном далее тексте программы в фигурных скобках, как правило, на
русском языке, за исключением тех случаев, когда в комментарии
указывается английскими буквами название процедуры в конце ее
описания или отмечается конец какого-либо блока, например,
end;{eigenvals}
end;{case}
end;{i}…
Замечание 1.4. Далее, при записи текстов программ все служебные слова будем выделять жирным шрифтом.
Текст программы lr1 для варианта 0
(t0 = 0, tt = 1, dt1 = 0.1, dt2 = 0.01)
1) program lr1; {Задача Коши: dx/dt =Ax, x(t0)=x0}
2) uses CRT, Graph;
3) type matr = array [1..2,1..2] of real; vect = array [1..2] of real;
4) const t0 = 0.0; tt = 1.0; dt1 = 0.1; dt2 = 0.01;
5) xmin = – 3; xmax = 3; ymin = – 3; ymax = 3;
6) grPath = ‘ d : \ turbo’;
7) e : matr = ((1,0),(0,1)); a : matr = ((-2,1),(0, -3));
8) x0 : vect = (2,2);
9) var rr, i, j, dv, mv, xc, yc : integer ;
10) aa, bb, t, kx, ky : real; ll, x, xx : vect;
11) function det(a : matr) : real;
12) begin
12
13)
14)
15)
16)
17)
18)
19)
20)
21)
22)
23)
24)
25)
26)
27)
28)
29)
30)
31)
32)
33)
34)
35)
36)
37)
38)
39)
40)
41)
42)
43)
44)
45)
46)
47)
48)
49)
50)
51)
52)
53)
54)
55)
56)
57)
58)
59)
det : = a[1,1]*a[2,2] – a[1,2]*a[2,1]
end;
procedure eigenvals(a : matr ; var ll : vect;
var aa, bb: real; var rr : integer);
{Нахождение собственных значений матрицы A};
var p, q, d : real;
begin
p : = – ( a[1,1] + a[2,2] ); q : = det(a);
d : = sqr(p)/4 – q; aa : = – p/2;
if d > 1E–6 then
begin
bb : = sqrt(d);
ll[1] : = aa + bb; ll[2] : = aa – bb;
rr : = 1
end
else if abs(d) < 1E–6 then
begin
ll[1] : = aa; ll[2] : = aa;
rr : = 2
end
else
begin
bb : = sqrt(–d);
rr : = 3
end
end; { eigenvals}
procedure cmatr(c : real; m: matr ; var cm : matr);
{Умножение матрицы на число}
var i, j : integer ;
begin
for i : = 1 to 2 do
for j : = 1 to 2 do
cm[i,j] : = c*m[i,j]
end;
procedure smatr(m1, m2 : matr ; var m : matr);
{Сложение матриц}
var i, j : integer ;
begin
for i : = 1 to 2 do
for j : = 1 to 2 do
m[i,j] : = m1[i,j] + m2[i,j]
end;
procedure pmatr(m: matr ; v : vect ; var p : vect);
{Умножение квадратной матрицы на матрицу-столбец}
var i, j : integer ; s : real ;
begin
for i : = 1 to 2 do
13
60)
begin
61)
s : = 0;
62)
for j : = 1 to 2 do
63)
s : = s + m[i,j]*v[j] ;
64)
p[i] : = s
65)
end
66) end;
67) procedure Coshy(a: matr ; x0: vect ; rr : integer ;
68) t0, t : real ; var x : vect);
69) {Решение задачи Коши; результат x – вектор составляющих
70) частного решения в момент t}
71) var ll1e, ll2e, aae, q1, q2, e1, e2 ,e3, expa : matr;
72) begin
73)
case rr of
74)
1: begin
75)
cmatr(–ll[1], e, ll1e); cmatr(–ll[2], e, ll2e);
76)
smatr(a, ll1e, q1); smatr(a, ll2e, q2);
77)
cmatr(exp(ll[1]*(t–t0))/(ll[1] – ll[2]),q2, e1);
78)
cmatr(exp(ll[2]*(t–t0))/(ll[2] – ll[1]),q1, e2);
79)
smatr(e1, e2, expa)
80)
end;
81)
2: begin
82)
cmatr(–ll[1], e, ll1e); smatr(a, ll1e, q1);
83)
cmatr (t–t0, q1, e1); smatr(e, e1, e2);
84)
cmatr(exp(ll[1]*( t–t0)), e2, expa)
85)
end;
86)
3: begin
87)
cmatr(cos(bb*( t–t0)), e, e1);
88)
cmatr(–aa, e, aae); smatr(a, aae, q1);
89)
cmatr(sin(bb*( t–t0))/bb, q1, e2);
90)
smatr(e1, e2, e3);
91)
cmatr(exp(aa*( t–t0)), e3, expa)
92)
end
93)
end; {case}
94)
pmatr(expa, x0, x)
95) end; { Coshy }
96) procedure CoshyResult(t: real; var xx: vect);
97) {Процедура введена для контроля вычислений}
98) begin
99)
xx[1] : = 4*exp(-2*t)-2*exp(-3*t);
100)
xx[2] : =2*exp(-3*t)
101) end;
102) {Далее следуют описания процедур и функций,
103) связанных с использованием графического режима}
104) function xscr (x: real): integer ;
105) {Переход от реальной координаты x к экранной}
106) begin
14
107)
108)
109)
110)
111)
112)
113)
114)
115)
116)
117)
118)
119)
120)
121)
122)
123)
124)
125)
126)
127)
128)
129)
130)
131)
132)
133)
134)
135)
136)
137)
138)
139)
140)
141)
142)
143)
144)
145)
146)
147)
148)
149)
150)
151)
152)
153)
xscr : = round( xc + x*kx)
end;
function yscr (y: real): integer ;
{Переход от реальной (математической) координаты y к экранной}
begin
yscr : = round( yc– y*ky)
end;
procedure myMoveTo(x,y: real);
{Перемещение графического указателя
в точку с реальными координатами x, y}
begin
MoveTo(xscr (x), yscr (y))
end;
procedure myLineTo(x,y: real);
{Построение отрезка из точки, в которой находится
графический указатель, в точку с реальными координатами x, y}
begin
LineTo(xscr (x), yscr (y))
end;
procedure OutStr (x: real);
{Вывод на экран значения переменной x в формате x : 4: 1}
var st : string ;
begin
Str(x:4:1, st); OutText(st)
end;
procedure osy;
{Построение осей координат и вывод значений
xmin, xmax, ymin, ymax около осей}
begin
kx : = (GetMaxX – 70)/(xmax – xmin);
ky : = (GetMaxY – 70)/(ymax – ymin);
xc : = round(abs(xmin*kx)) + 30;
yc : = round(abs(ymax*ky)) + 30;
myMoveTo(xmin, 0); myLineTo(xmax,0);
OutTextXY(GetX, GetY – 3, ‘ >x1’);
myMoveTo(0, ymin); myLineTo(0, ymax);
OutTextXY(GetX – 3, GetY, ‘ ^x2’);
myMoveTo(xmin , – 0.1); OutStr(xmin);
myMoveTo(xmax – 0.2, – 0.1); OutStr(xmax);
myMoveTo( – 0.4, ymin+0.1); OutStr(ymin);
myMoveTo( – 0.4, ymax-0.1); OutStr(ymax)
end; {osy}
begin {Основная программа}
ClrScr ;
{Вычисления}
writeln(‘ Матрица A
Вектор x0’);
for i : = 1 to 2 do
15
154)
begin
155)
for j : = 1 to 2 do write (a[i,j] : 10 : 4);
156)
write(x0[i] : 15 : 4);
157)
writeln
158)
end;
159)
eigenvals(a, ll, aa, bb, rr);
160)
writeln(‘rr=’, rr);
161)
if rr < 3 then
162)
writeln(‘ ll[1] =’, ll[1] : 8 : 4, ‘ll[2] =’, ll[2] : 8 : 4)
163)
else
164)
writeln(‘ ll[1] =’, aa : 8 : 4,‘ +’, bb : 8 : 4, ‘i ll[2] =’,
165)
aa : 8: 4,‘ – ’, bb : 8 : 4, ‘i’);
166)
t : = t0; CoshyResult(t, xx);
167)
writeln(‘ t
x1
x2
xx1
xx2 ’);
168)
writeln(t: 8:2, x0[1]:12:4, x0[2]:12:4, xx[1]:12:4, xx[2]:12:4);
169)
while t <= tt do
170)
begin
171)
t : =t + dt1;
172)
Coshy(a, x0, rr, t0, t, x);
173)
CoshyResult(t, xx);
174)
writeln(t:8:2, x[1]:12:4, x[2]:12:4, xx[1]:12:4, xx[2]:12:4)
175)
end;
176)
readln;
177) {График}
178)
dv : = detect;
179)
InitGraph(dv, mv, grPath);
180)
if GraphResult < > 0 then Halt(1);
181)
SetBkColor(1);
182)
osy;
183)
SetColor(5);
184)
t : = t0;
185)
myMoveTo(x0[1], x0[2]);
186)
while t < = tt do
187)
begin
188)
t : = t + dt2;
189)
Coshy(a, x0, rr, t0, t, x);
190)
myLineTo(x[1], x[2])
191)
end;
192)
readln;
193)
CloseGraph
194) end.
Результаты решения
I. Результаты, связанные с вычислениями:
Матрица А
Вектор х0
-2.0000 1.0000
2.0000
0.0000 -3.0000
2.0000
rr  1
ll 1  3.0000
ll 2  2.0000
t
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
x1
2.0000
1.7933
1.5837
1.3821
1.1949
1.0253
0.8742
0.7415
0.6262
0.5268
0.4418
x2
2.0000
1.4816
1.0976
0.8131
0.6024
0.4463
0.3306
0.2449
0.1814
0.1344
0.0996
xx1
2.0000
1.7933
1.5837
1.3821
1.1949
1.0253
0.8742
0.7415
0.6262
0.5268
0.4418
xx2
2.0000
1.4816
1.0976
0.8131
0.6024
0.4463
0.3306
0.2449
0.1814
0.1344
0.0996
II. Результаты, связанные с построением графика (рис.1.1):
Рис. 1.1
Замечание 1.5. Сведения о стандартных подпрограммах, связанных с работой в графических режимах, содержатся, например,
в книгах [7, 8].
После выполнения лабораторной работы студенты сдают отчет. Отчет должен содержать следующие разделы:
16
17
Окончание таблицы
1. Постановка задачи.
2. Текст программы в системе MathCad.
3. Результаты решения на языке Паскаль.
№
варианта
12
Система дифференциальных
уравнений
13
 x1   x1

 x2  3x1  x2
Начальные условия
14
x1 0   2
x 2 0   2
 x1  2 x1  x2

 x 2  x1  2 x2
x1 0   2
x 2 0   2
15
x1 0   2
x 2 0   2
 x1   x1

 x 2  x1  x2
x1 0   2
x 2 0   2
16
x1 0   2
x 2 0   2
 x1  2 x1  x2

 x 2  4 x1  2 x2
x1 0   2
x 2 0   1
17
 x1  2 x1  x2

 x 2   x2
x1 0   2
x 2 0   2
Варианты для выполнения лабораторной работы
t  0,1, t  0.1
 x1  5 x1  x2

 x2  2 x2
Начальные условия
x1 0   2
x 2 0   2
x1 0   2
x 2 0   2
№
варианта
1
Система дифференциальных
уравнений
2
 x1  2 x1

 x 2  x1  5 x2
3
 x1   x1  3x2

 x 2   x2
4
 x1  2 x1  x2

 x 2   x1  2 x2
x1 0   2
x 2 0   2
18
5
 x1   x1  x2

 x2   x2
x1 0   2
x 2 0   2
 x1   x1  2 x2

 x 2   x1  3x2
19
6
 x1  2 x1  4 x2

 x 2  x1  2 x2
x1 0   1
x 2 0   2
 x1  3x1  x2

 x 2  x1  x2
x1 0    1
x 2 0   1
20
7
 x1   x1

 x 2  x1  2x2
x1 0   2
x 2 0   2
 x1  2 x1

 x 2  x1  x2
x1 0   2
x 2 0   2
8
 x1  3x1  x2

 x 2  2 x1  x2
x1 0   2
x 2 0   2
9
 x1   x1  x2

 x 2   x1  3x2
x1 0   1
x 2 0   1
10
 x1   x1  x2

 x 2  2x2
x1 0   2
x 2 0   2
11
 x1  2 x1

 x 2  2 x1  x2
 x1   x1  2 x2

 x 2  2 x2
18
x1 0   2
x 2 0   2
19
x1 0   2
x 2 0   2
ЛАБОРАТОРНАЯ РАБОТА № 2
Решение задачи Коши для простейшей системы двух линейных
неоднородных дифференциальных уравнений 1-го порядка
с постоянными коэффициентами
Необходимые теоретические сведения
В лабораторной работе № 1 было рассмотрено решение задачи
Коши для системы, векторная запись которой имеет вид
x  Ax .
В данной лабораторной работе будем рассматривать решение
задачи Коши для несколько более сложной системы вида
x  Ax  v .
Как будет показано далее, эта задача легко может быть сведена
к предыдущей.
Итак, в работе рассматривается система дифференциальных
уравнений
 x1  a11 x1  a12 x2  v1,
(2.1)

 x2  a21 x1  a22 x2  v2 .
Требуется найти частное решение этой системы, удовлетворяющее начальным условиям
(2.2)
x1 0   x10 , x2 0   x20 ,
т. е. найти решение задачи Коши.
Задача (2.1), (2.2) в векторной форме имеет вид:
найти частное решение векторного уравнения
(2.3)
x  Ax  v ,
удовлетворяющее начальному условию
(2.4)
x0   x 0 .
Для сведения этой задачи к задаче, решавшейся в лабораторной
работе № 1, преобразуем правую часть векторного уравнения (2.3):
Ax  v  Ax  AA1v  A x  A1v  A x  c   Ay .
Введем обозначения
c  A1v,
y  x  c.
Ясно, что y  x  c  x  0  x , левая часть векторного уравнения (2.3) равна y , а правая Ay и система (2.3) принимает вид
y  Ay.
(2.5)

20

Начальные условия
y ( 0)  x ( 0)  c  x 0  c .
(2.6)
Получившаяся задача (2.5), (2.6) имеет вид задачи (1.9), (1.10),
рассмотренной в лабораторной работе № 1. Частное решение уравнения (2.5) находится по формуле
y t   e At y 0 ,
следовательно
xt   c  e At x 0  c ,
откуда
(2.7)
xt   e At x 0  c  c,
где
c  A1v.
(2.8)
1
Как нетрудно проверить, матрица A имеет вид
1  a22  a12 

.
(2.9)
A 1 
det A   a21 a11 




Выполнение лабораторной работы № 2
Вариант 0
1. Постановка задачи
Дана система дифференциальных уравнений
 x1  2 x1  x2  1,

 x2  3x2  1.
Требуется найти частное решение системы, удовлетворяющее
начальным условиям
x1 0   2, x 2 0   2 .
2. Решение в системе MathCad
Для решения задачи нам потребуются следующие данные:
 2 1
а) A  
 – матрица коэффициентов правой части системы;
 0  3
1
б) v    – постоянный вектор в правой части системы;
  1
 2
в) x0    – вектор начальных условий;
 2
г) верхняя граница промежутка интегрирования T , например,
T  1;
21
д) формулы составляющих искомого частного решения;
1
1
xx1 t   4e  2t  2 e 3t  ;
3
3
1
1
xx2 t   2 e  3t  .
3
3
Эти формулы используются для проверки правильности работы программы.
Текст программы
ORIGIN : 1
Единичная матрица
1 0

E : 
0 1
Матрица коэффициентов правой части системы
 2 1 

A : 
 0  3
Столбец свободных членов
 1
v :  
  1
Начальные условия
 2
x0 :  
 2
Верхняя граница промежутка интегрирования T : 1
Шаг аргумента h : 0.1
T
Число шагов N : round  
N  10
h
Собственные значения матрицы А
  2
λ : eigenvalsA 
λ   
  3
 : Re1 
 : Im1 
  2
0
Нахождение матричной экспоненты:
λ2t
e 1 t
EXP 1t  :
A  λ 2 E   e
A  λ1E 
1   2
 2  1


A  αE
sin β t 
EXP 3t  : e α t  Ecos β t  
β


EXPAt  : if λ1  λ 2 , EXP 2t , if β  0, EXP1t , EXP 3t 
Нахождение вектора-столбца c:
  0.333
c : A 1  v
c  

 0.333 
Начальные значения искомых функций:
X11 : x 01
X 21 : x 0 2
i : 1.. N 1
Значение аргумента
t i : h i  1
Частное решение:
 X1i 
 : EXPAt i x 0  c  c

 X 2i 
Решение, полученное аналитически:
7
1
XX1i : 4e  2 t i  e 3t i 
3
3
7 3t i 1
XX2 i : e
 .
3
3
ti =
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
X1i =
2
1.88
1.734
1.58
1.428
1.284
1.152
1.034
0.929
0.838
0.759
X2i =
2
1.395
0.947
0.615
0.369
0.187
0.052
-0.048
-0.122
-0.177
-0.217
EXP 2t  : e λ1t E  A  λ1E  t 
22
23
XX1i =
2
1.88
1.734
1.58
1.428
1.284
1.152
1.034
0.929
0.838
0.759
XX2i =
2
1.395
0.947
0.615
0.369
0.187
0.052
-0.048
-0.122
-0.177
-0.217
3. Решение на языке Паскаль
Для выполнения работы на языке Паскаль используется приводимая далее программа lr2, предназначенная для решения задачи
Коши для простейшей системы двух линейных неоднородных уравнений с постоянными коэффициентами.
Исходные данные
Исходные данные задаются в разделе описаний констант
(см. далее текст программы, строки 4–8):
t 0 – нижняя граница промежутка интегрирования;
tt – верхняя граница промежутка интегрирования;
dt1 – шаг аргумента t при выводе результатов решения;
dt 2 – шаг аргумента t при построении графика;
x min, x max – нижняя и верхняя границы изменения переменной;
y min, y max – нижняя и верхняя границы изменения переменной x2 ;
a – двумерный массив элементов матрицы А (матрицы коэффициентов при неизвестных в правой части системы);
x 0 – одномерный массив, вектор-столбец начальных условий;
grPath – константа строкового типа, путь к библиотеке графических подпрограмм (3-й параметр в обращении к стандартной процедуре InitGraph, в приведенном тексте grPath = ‘d : \ turbo’);
процедура CoshyResult введена для проверки правильности работы программы, задает точные формулы составляющих x1 (t ) ,
x2 (t ) частного решения, полученные заранее; формулы для x1 (t ) ,
x2 (t ) задаются с помощью одномерного массива xx с составляющими xx1, xx2 соответственно (см. далее текст программы, строки 116–121).
Результаты решения
I. Результаты, связанные с вычислениями
1. Исходные данные: матрица А (массив А), вектор v (массив
v) и начальный вектор x 0 (массив x0 ), см. строки 172–178.
2. Число rr и собственные значения матрицы А, см. строки
180–185:
24
а) число rr :
если собственные значения матрицы А вещественны и различны, то rr  1 ;
если собственные значения матрицы А совпадают, то rr  2 ;
если собственные значения матрицы А комплексно-сопряженные, то rr  3 ;
б) собственные значения матрицы А:
если собственные значения λ1 и λ 2 матрицы А вещественны,
то выводится результат
ll 1  λ1 ll 2  λ 2 ;
если собственные значения λ1 и λ 2 матрицы А комплексносопряженные λ 1,2  α  βi, то выводится результат
ll 1  α  βi ll 2  α  βi;
3) таблица результатов, состоящая из пяти столбцов, см. строки 187–188, 194:
 в 1-м столбце выводятся значения переменной t от t 0 до tt
с шагом dt1;
 во 2-м и 3-м столбцах выводятся соответствующие значения
составляющих x1 и x2 массива x – значения составляющих x1 t 
и x2 t  частного решения (они получаются в результате выполнения
процедуры Coshy);
 в 4-м и 5-м столбцах выводятся для контроля значения составляющих xx1 и xx2 массива xx – значения составляющих частного решения, формулы для которых задаются в процедуре
CoshyResult.
II. Результаты, связанные с построением графика
Строятся оси координат, на которых указываются значения
x min, x max, y min, y max, см. строки 202, 152–168.
Строится траектория движения объекта из точки с координатами x10 , x20 по закону x1 t , x2 t  путем построения отрезков прямых, соединяющих точки с координатами x1, x2, построенные
последовательно в моменты времени из промежутка t 0; tt  с шагом
dt 2 , см. строки 205–211.
25
Замечание 2.1. Программа lr2 составлена на основе программы lr1. Для удобства строки, совпадающие со строками программы
lr1, помечены знаком «*» (см. далее текст программы).
Текст программы lr2 для варианта 0
(t0=0, tt=1, dt1=0.1, dt2=0.01)
1) program lr2; {Задача Коши: dx/dt =Ax+v, x(t0)=x0}
2) uses CRT, Graph;
3) type matr = array [1..2,1..2] of real; vect = array [1..2] of real;
4) const t0 = 0.0; tt =1.0; dt1 = 0.1; dt2 =0.01;
5) xmin = – 0.5; xmax = 2.5; ymin = – 1; ymax = 2;
6) grPath = ‘d : \ turbo’;
7) e : matr = ((1,0),(0,1)); a : matr = ((–1,2),(0, –2));
8) v: vect = (1, –1); x0 : vect = (2,2);
9) var rr, i, j, dv, mv, xc, yc : integer ;
10) aa, bb, t, kx, ky : real; ll, x, xx : vect;
11) function det(a : matr) : real;
12) begin
13) det : = a[1,1]*a[2,2] – a[1,2]*a[2,1]
14) end;
15) procedure eigenvals(a : matr ; var ll : vect;
16) var aa, bb: real; var rr : integer);
17) {Нахождение собственных значений матрицы A}
18) var p, q, d : real;
19) begin
20)
p : = – ( a[1,1] + a[2,2] ); q : = det(a);
21)
d : = sqr(p)/4 – q; aa : = –p/2;
22)
if d > 1E–6 then
23)
begin
24)
bb : = sqrt(d);
25)
ll[1] : = aa + bb ; ll[2] : = aa – bb;
26)
rr : = 1
27)
end
28)
else if abs(d) < 1E–6 then
29)
begin
30)
ll[1] : = aa; ll[2] : = aa;
31)
rr : = 2
32)
end
33)
else
34)
begin
35)
bb : = sqrt(–d);
36)
rr := 3
37)
end
38) end; { eigenvals}
26
*
*
……*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
….. *
*
*
*
*
*
*
*
*
*
*
*
*
*
39) procedure cmatr(c : real; m: matr ; var cm : matr);
40) {Умножение матрицы на число}
41) var i, j : integer ;
42) begin
43)
for i: = 1 to 2 do
44) for j : = 1 to 2 do
45)
cm[i,j] : = c*m[i,j]
46) end;
47) procedure smatr(m1, m2 : matr ; var m : matr);
48) {Сложение матриц}
49) var i, j : integer ;
50) begin
51)
for i : = 1 to 2 do
52)
for j : = 1 to 2 do
53)
m[i,j] : = m1[i,j] + m2[i,j]
54) end;
55) procedure pmatr(m: matr ; v : vect ; var p : vect);
56){Умножение квадратной матрицы на матрицу-столбец}
57) var i, j : integer ; s : real ;
58) begin
59) for i : = 1 to 2 do
60)
begin
61)
s : = 0;
62)
for j : = 1 to 2 do
63)
s : = s + m[i,j]*v[j] ;
64)
p[i] : = s
65)
end
66) end;
67) procedure cvect(c : real; v: vect ; var cv : vect);
68) {Умножение вектора на число}
69) var i : integer ;
70) begin
71)
for i : = 1 to 2 do cv[i] : = c*v[i]
72) end;
73) procedure svect(v1, v2 : vect ; var v : vect);
74) {Сложение векторов}
75) var i : integer ;
76) begin
77)
for i : = 1 to 2 do v[i] : = v1[i] + v2[i]
78) end;
79) procedure Coshy(a: matr ; v, x0: vect ; rr : integer ;
80) t0, t : real ; var x : vect);
81) {Решение задачи Коши; результат x – вектор составляющих
82) частного решения в момент t}
83) var ll1e, ll2e, aae, q1, q2, e1, e2 ,e3, expa,
27
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
84) a1, a2 : matr; c, c1, x0c,xc : vect;
85) begin
*
86)
a1[1,1] : = a[2,2]; a1[1,2] : = –a[1,2];
87)
a1[2,1] : = –a[2,1]; a1[2,2] : = a[1,1];
88)
cmatr(1/det(a), a1, a2);
89)
pmatr(a2, v, c);
90)
svect(x0, c, x0c);
91) case rr of
*
92) 1: begin
*
93)
cmatr(–ll[1], e, ll1e); cmatr(–ll[2], e, ll2e);
*
94)
smatr(a, ll1e, q1); smatr(a, ll2e, q2);
*
95)
cmatr(exp(ll[1]*(t–t0))/(ll[1] – ll[2]),q2, e1);
*
96)
cmatr(exp(ll[2]*(t–t0))/(ll[2] – ll[1]),q1, e2);
*
97)
smatr(e1, e2, expa)
*
98)
end;
*
99) 2: begin
*
100)
cmatr(–ll[1], e, ll1e); smatr(a, ll1e, q1);
*
101)
cmatr (t–t0, q1, e1); smatr(e, e1, e2);
*
102)
cmatr(exp(ll[1]*( t–t0)), e2, expa)
*
103)
end;
*
104) 3: begin
*
105)
cmatr(cos(bb*( t–t0)), e, e1);
*
106)
cmatr(–aa, e, aae); smatr(a, aae, q1);
*
107)
cmatr(sin(bb*( t – t0))/bb, q1, e2);
*
108)
smatr(e1, e2, e3);
*
109)
cmatr(exp(aa*( t – t0)), e3, expa)
*
110)
end
*
111) end; {case}
*
112) pmatr(expa, x0c, xc);
113) cvect( –1, c, c1);
114) svect(xc, c1, x)
115) end; { Coshy }
*
116) procedure CoshyResult(t: real; var xx: vect);
*
117) {Процедура введена для контроля вычислений}
*
118) begin
*
119) xx[1] : = 7*exp(– t) – 5*exp(– 2*t);
120) xx[2] : =2.5*exp(– 2*t) – 0.5
121) end;
*
122) {Далее следуют описания процедур и функций,
*
123) связанных с использованием графического режима}
*
124) function xscr (x: real): integer ;
*
125){Переход от реальной (математической) координаты x к экранной} *
126) begin
*
127) xscr : = round( xc + x*kx)
*
128) end;
*
28
129) function yscr (y: real): integer ;
*
130){Переход от реальной (математической) координаты y к экранной} *
131) begin
*
132)
yscr : = round( yc– y*ky)
*
133) end;
*
134) procedure myMoveTo(x,y: real);
*
135) {Перемещение графического указателя
*
136) в точку с реальными координатами x, y}
*
137) begin
*
138) MoveTo(xscr (x), yscr (y))
*
139) end;
*
140) procedure myLineTo(x,y: real);
*
141) {Построение отрезка из точки, в которой находится
*
142) графический указатель, в точку с реальными координатами x, y} *
143) begin
*
144) LineTo(xscr (x), yscr (y))
*
145) end;
*
146) procedure OutStr (x: real);
*
147) {Вывод на экран значения переменной x в формате x : 4: 1}
*
148) var st : string ;
*
149) begin
*
150) Str(x:4:1, st); OutText(st)
*
151) end;
*
152) procedure osy;
*
153) {Построение осей координат и вывод значений
*
154) xmin, xmax, ymin, ymax около осей}
*
155) begin
*
156) kx : = (GetMaxX – 70)/(xmax – xmin);
*
157) ky : = (GetMaxY – 70)/(ymax – ymin);
*
158) xc : = round(abs(xmin*kx)) + 30;
*
159) yc : = round(abs(ymax*ky)) + 30;
*
160) myMoveTo(xmin, 0); myLineTo(xmax,0);
*
161) OutTextXY(GetX, GetY – 3, ‘ >x1’);
*
162) myMoveTo(0, ymin); myLineTo(0, ymax);
*
163) OutTextXY(GetX – 3, GetY, ‘ ^x2’);
*
164) myMoveTo(xmin , 0.12); OutStr(xmin);
165) myMoveTo(xmax, 0.12); OutStr(xmax);
166) myMoveTo( – 0.25, ymin); OutStr(ymin);
167) myMoveTo( – 0.25, ymax); OutStr(ymax)
168) end; {osy}
*
169) begin {Основная программа}
*
170) ClrScr ;
*
171) {Вычисления}
*
172)
writeln(‘ Матрица A
Вектор v
Вектор x0’);
173)
for i : = 1 to 2 do
*
29
174)
begin
175)
for j : = 1 to 2 do write (a[i,j] : 10 : 4);
176)
write(v[i] : 15 : 4, x0[i] : 15 : 4);
177)
writeln
178)
end;
179)
eigenvals(a, ll, aa, bb, rr);
180)
writeln(‘rr=’, rr);
181)
if rr < 3 then
182)
writeln(‘ ll[1] =’, ll[1] : 8 : 4, ‘ ll[2] =’, ll[2] : 8 : 4)
183)
else
184)
writeln(‘ ll[1] =’, aa : 8 : 4,‘ +’, bb : 8 : 4, ‘i ll[2] =’,
185)
aa : 8: 4,‘ – ’, bb : 8 : 4, ‘i’);
186)
t : = t0; CoshyResult(t, xx);
187)
writeln(‘ t
x1
x2
xx1
xx2 ’);
188)
writeln(t: 8:2, x0[1]:12:4, x0[2]:12:4, xx[1]:12:4, xx[2]:12:4);
189)
while t <= tt do
190)
begin
191)
t : =t + dt1;
192)
Coshy(a, v, x0, rr, t0, t, x);
193)
CoshyResult(t, xx);
194)
writeln(t:8:2, x[1]:12:4, x[2]:12:4, xx[1]:12:4, xx[2]:12:4)
195)
end;
196)
readln;
197)
{График}
198)
dv : = detect;
199)
InitGraph(dv, mv, grPath);
200)
if GraphResult < > 0 then Halt(1);
201)
SetBkColor(1);
202)
osy;
203)
SetColor(5);
204)
t : = t0;
205)
myMoveTo(x0[1], x0[2]);
206)
while t < = tt do
207)
begin
208)
t : = t + dt2;
209)
Coshy(a, v, x0, rr, t0, t, x);
210)
myLineTo(x[1], x[2])
211)
end;
212)
readln;
213)
CloseGraph
214) end.
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
Результаты решения
I. Результаты, связанные с вычислениями
Матрица A
Вектор v
–1.0000
1.0000
1.0000
0.0000 –3.0000
–1.0000
rr  1
ll 1  3.0000
ll 2  2.0000
t
x1
x2
0.00
2.0000
2.0000
0.10
1.8797
1.3952
0.20
1.7341
0.9472
0.30
1.5799
0.6153
0.40
1.4279
0.3695
0.50
1.2842
0.1873
0.60
1.1524
0.0524
0.70
1.0340
– 0.0476
0.80
0.9292
– 0.1217
0.90
0.8377
– 0.1765
1.00
0.7585
– 0.2172
Вектор x0
2.0000
2.0000
xx1
2.0000
1.8797
1.7341
1.5799
1.4279
1.2842
1.1524
1.0340
0.9292
0.8377
0.7585
xx2
2.0000
1.3952
0.9472
0.6153
0.3695
0.1873
0.0524
– 0.0476
– 0.1217
– 0.1765
– 0.2172
II. Результаты, связанные с построением графика (рис. 2.1)
Рис. 2.1
После выполнения данной работы студенты сдают отчет. Отчет должен содержать следующие разделы:
30
31
Окончание таблицы
1. Постановка задачи.
2. Текст программы в системе MathCad.
3. Результаты решения на языке Паскаль.
№
варианта
12
Система дифференциальных
уравнений
13
 x1   x1  1

 x 2  3x1  x2  1
Начальные условия
14
x1 0   2
x 2 0   2
 x1  2 x1  x2  1

 x 2  x1  2 x2  1
x1 0   2
x 2 0   2
15
x1 0   2
x 2 0   2
 x1   x1  1

 x 2  x1  x2  1
x1 0   2
x 2 0   2
16
 x1   x1  3x2  1

 x 2   x2  1
x1 0   2
x 2 0   2
 x1  2 x1  x2  1

 x 2  4 x1  2 x2  1
x1 0   2
x 2 0   1
17
 x1  2 x1  x2  1

 x 2   x1  2 x2  1
x1 0   2
x 2 0   2
 x1  2 x1  x2  2

x 2   x2  1

x1 0   2
x 2 0   2
18
 x1   x1  x2  1

 x 2   x2  1
x1 0   2
x 2 0   2
 x1   x1  2 x2  1

 x 2   x1  3x2  1
19
 x1  2 x1  4 x2  1

 x 2  x1  2 x2  1
x1 0   1
x 2 0   2
 x1  3x1  x2  2

 x 2  x1  x2  1
x1 0    1
x 2 0   1
20
 x1  2 x1  2

 x 2  x1  x2  1
x1 0   2
x 2 0   2
Варианты для выполнения лабораторной работы
t  0,1, t  0.1
№
варианта
1
2
3
4
5
6
Система дифференциальных
уравнений
 x1   x1  2 x2  1

 x 2  2 x2  1
 x1  2 x1  1

 x 2  x1  5 x2  1
x1 0   2
x 2 0   2
7
 x1   x1  1

 x 2  x1  2 x2  2
8
 x1  3x1  x2  1

 x 2  2 x1  x2  1
9
 x1   x1  x2  1

 x 2   x1  3x2  2
10
 x1   x1  x2  1

 x 2  2 x2  2
x1 0   2
x 2 0   2
11
 x1  2 x1  1

 x 2  2 x1  x2  1
x1 0   2
x 2 0   2
32
 x1  5 x1  x2  1

 x 2  2 x2  1
x1 0   2
x 2 0   2
x1 0   1
x 2 0   1
33
Начальные условия
x1 0   2
x 2 0   2
x1 0   2
x 2 0   2
x1 0   2
x 2 0   2
ЛАБОРАТОРНАЯ РАБОТА № 3
Линейная задача оптимального быстродействия. Построение
траекторий для системы двух дифференциальных уравнений
1-го порядка со скалярным управлением u
Теоретические сведения, необходимые для выполнения
лабораторных работ № 3, 4
1. Постановка линейной задачи оптимального быстродействия
При выполнении данной лабораторной работы рассматривается управляемый объект с числом фазовых координат n  2 и числом
управляющих параметров r  1. Система уравнений движения этого
объекта имеет вид:
 x1  a11 x1  a12 x 2  b1u ,
(3.1)

 x 2  a 21 x1  a 22 x 2  b 2 u .
Здесь a11 , a12 , a21 , a22 , b1 , b2 – заданные постоянные числа.
Требуется найти управление u  u (t ) , удовлетворяющее ограничению
u 1
(3.2)
и переводящее объект из заданного начального состояния
x1 0   x10 , x2 0   x20
(3.3)
в начало координат
x1 T   0, x2 T   0
(3.4)
за наименьшее возможное время T . Это – задача об оптимальном
быстродействии.
в начало координат
x(T )  0
за наименьшее возможное время T .
(3.8)
3. Допустимое управление
Допустимым управлением называется любая кусочнонепрерывная функция u  u (t ) , t  t0 ;t1  со значениями, удовлетворяющими ограничению (3.2) (или (3.6)), непрерывная справа в точках разрыва и непрерывная на концах отрезка t0 ;t1 , на котором она
задана.
4. Сопряженные переменные. Сопряженная система дифференциальных уравнений
Наряду с переменными x1 , x2 (фазовыми координатами объекта) будем рассматривать вспомогательные переменные 1 ,  2 (координаты вектора  ), которые называются сопряженными переменными. Закон их изменения задается векторным дифференциальным уравнением
   AT 

(3.9)
или системой дифференциальных уравнений
  1  a111  a21 2 ,
(3.10)

 2  a12 1  a22  2 .

Система (3.10) называется сопряженной системой по отношению к исходной системе (3.1).
5. Принцип максимума Л. С. Потряхина для линейной задачи
оптимального быстродействия
Теорема 3.1
Пусть u (t ) , 0  t  T – допустимое управление, переводящее
2. Постановка линейной задачи оптимального быстродействия в векторной форме
Пусть движение управляемого объекта описывается векторным уравнением
(3.5)
x  Ax  bu .
Требуется найти управление, удовлетворяющее ограничению
u 1
(3.6)
и переводящее объект из заданного начального состояния
x(0)  x 0
(3.7)
ветствующая траектория (см. (3.5)), так что x(0)  x 0 , x(T )  0 . Для
оптимальности по быстродействию управления u (t ) и траектории
x(t ) необходимо существование такой ненулевой непрерывной вектор-функции (t )  (1 (t ),  2 (t )) , соответствующей функциям u (t )
и x(t ) (см. (3.9)), что для любого момента t 0  t  T выполняется
условие
34
35
фазовую точку из положения x 0 в начало координат, а x(t ) – соот-
u (t )  sign(b1 1 (t )  b2  2 (t )) .
(3.10)
Следствие 3.1
Оптимальное управление в линейной задаче оптимального быстродействия имеет релейный характер, т. е. принимает два значения: u  1, u  1.
6. Теорема А. А. Фельдбаума о числе переключений
Напомним, что собственными значениями матрицы А называa 
a12
 0.
ются корни уравнения det(A  E )  0 или 11
a 21
a 22  
Теорема 3.2
Рассматривается управляемый объект (3.5). Пусть собственные
значения 1 ,  2 матрицы А (см. (3.5)) вещественны. Тогда оптимальное по быстродействию управление имеет не более одного момента переключения, т. е. не более двух интервалов постоянства.
Следствие 3.2
Оптимальное по быстродействию управление в случае вещественных значений 1 и  2 матрицы А имеет график одного из видов,
показанных на рис. 3.1.
7. Принцип максимума – необходимое и достаточное условие
оптимальности
Теорема 3.3
Пусть в системе (3.5) выполнено условие:
(3.12)
rank b, Ab  2 .
Тогда для оптимальности управления u (t ) необходимо и достаточно, чтобы оно удовлетворяло принципу максимума (см. теорему (3.1).
Условие (3.12) называется условием общности положения. Для
r > 1 оно имеет более сложный вид.
Следствие 3.3
Пусть найдено управление вида (3.11), переводящее объект
(3.5) из состояния x(0)  x 0 в начало координат x(T )  0 и выполнено условие общности положения (3.12). Тогда u (t ) – управление,
оптимальное в смысле быстродействия.
8. План решения линейной задачи оптимального быстродействия [2]
Задача I. Найти решение (t ) сопряженной системы (3.9) при
произвольном начальном векторе (0)   0 .
Задача II. Зная нетривиальное решение системы (3.9), найти
управление u (t ) , удовлетворяющее условию максимума (3.11).
Задача III. Зная управление u (t ) , найти соответствующую тра-
Рис. 3.1
екторию, исходящую из заданной начальной точки x 0 (найти решение системы (3.5)).
Известно, что задачи I, II, III можно решить и они имеют единственное решение. Таким образом, при разных начальных векторах
 0 будут получаться разные траектории, исходящие из x 0 . Если
удастся найти такой начальный вектор  , при котором траектория
проходит через начало координат, то эта траектория и соответствующее управление будут оптимальными.
Задача IV. Найти начальное значение  0 , при котором соответствующая траектория приходит в начало координат.
Как указано в книге [2], точное решение задачи IV неизвестно
и вряд ли возможно, поэтому вектор  0 находят приближенно.
36
37
Замечание 3.1. Решение задач I, II, III рассмотрено в лабораторной работе № 3. Приближенное решение задачи IV осуществляется в лабораторной работе № 4.
9. Теоремы существования и единственности для линейной задачи оптимального быстродействия
Теорема 3.4 (теорема единственности)
Пусть u1 (t ) и u 2 (t ) – два оптимальных управления, заданные
соответственно на отрезках 0  t  T1 и 0  t  T2 и переводящие
точку x 0 в начало координат. Тогда эти управления совпадают, т. е.
T1  T2 и u1 (t )  u 2 (t ) на 0; T1 .
Областью управляемости объекта (3.5) называется множество
всех точек x 0 фазовой плоскости X, из которых возможно при помощи какого-либо допустимого управления попасть в начало координат.
Теорема 3.5 (теорема существования)
Для любой точки x0 , принадлежащей области управляемости,
существует оптимальное по быстродействию управление, переводящее точку x 0 в начало координат (если существует допустимое
управление, то существует и оптимальное управление).
Теорема 3.6
Если выполнено условие (3.12) и матрица А устойчива, то есть
Re  i  0 (i  1,2) , то область управляемости совпадает с фазовой
плоскостью X.
Выполнение лабораторной работы № 3
Вариант 0
1. Постановка задачи и проверка условия общности положения
Пусть движение управляемого объекта описывается системой
дифференциальных уравнений
 x1  2 x1  x2  u
.

 x 2  3x2  u
На управляющий параметр u наложено ограничение
38
u  1.
Начальное состояние объекта
x1 (0)  2 , x 2 (0)  2 .
Требуется построить N различных траекторий, исходящих из
заданного начального состояния при t  0; TT  (т. е. найти решение
задач I, II, III, сформулированных в п. 8). Рассмотреть, например,
случаи:
TT  1;
N  3,
TT  1;
N  5,
TT  1;
N  10,
TT  2 ;(*)
N  100,
TT  2 ;
N  1000,
TT  3 .
N  1000,
Подобрать N и TT так, чтобы множество точек фазовой плоскости, образующих эти N траекторий, содержало начало координат
(в лабораторной работе № 4 из этих N траекторий будет выбрана
траектория, для которой расстояние от некоторой точки траектории
до начала координат будет наименьшим – это приближенное решение задачи IV, п. 8).
Проверка условия общности положения (3.12)
  2 1
 1
 1  3
A  
b, Ab  
 ,   3  3  6  0 ,
 , b    ,
  3 0
  1
  1  3
условие (3.12) выполнено.
2. Решение задачи
Сущность метода решения задачи
Нужно построить N различных траекторий, исходящих из заданной начальной точки. Для этого сначала нужно построить N
различных начальных сопряженных векторов  0 .
1. Выбор N различных начальных сопряженных векторов Ψ0
для построения N различных траекторий
Замечание 3.2. Сопряженная система (3.10) является однородной, следовательно, если (1 (t ),  2 (t )) – решение (3.10), то
(k1 (t ), k 2 (t )) – также решение (3.10), где k  const. Если k  0 , то
управление (3.11) для этих двух решений сопряженной системы будет одним и тем же, т. е. длина вектора (t )  (1 (t ),  2 (t )) не име39
ет значения, а важно только направление. Можно считать для определенности, что  0  (t 0 )  1 .
Нужно выбрать N различных сопряженных начальных векторов  0 . Обозначим
2
 .
(3.13)
N
Тогда в качестве начальных сопряженных векторов  0 можно
брать векторы с координатами
(cos l , sinl ),
2
где l  1, 2, ..., N . Например, для N  3 получается   , для N  8
3
2 
получается  
 .
8 4
2. Построение N различных траекторий для l  1, 2, ..., N ,
t  0, TT  , где TT – заданное число (далее считаем t0  0 ).
Таким образом, строится N различных траекторий, исходящих
из начальной точки (3.3), t  0; TT . Если при рассматриваемых N
и TT множество точек, образующих эти траектории, еще не содержит
начало координат – точку 0, то N и TT нужно увеличить, рассмотрев
случаи, приведенные в (*) или аналогичные им. Если множество точек, образующих траектории, содержит точку 0, то решение задачи
закончено (траектории «показывают» начало координат).
Решение задачи на языке Паскаль
Для выполнения данной работы на языке Паскаль используется приводимая далее программа lr3, предназначенная для построения N различных траекторий системы (3.1), исходящей из начальной точки (3.3) при t  0; TT , где TT – заданное число.
Исходные данные
Пусть l – номер сопряженного вектора  0  (0) и, следовательно, l – номер траектории, соответствующей этому сопряженному вектору. Для каждого l  1, 2, ..., N строится начальный сопряженный вектор
 (0)  (cos l , sinl),
(3.14)
где  находится по формуле (3.13).
Находится решение (t ) сопряженной системы (3.9) с начальным условием (3.14) на 0; TT , как это делалось в лабораторной работе № 1 (п. 8, задача I).
Находится управление u(t ) , соответствующее решению (t ) ,
удовлетворяющее условию (3.11) (п. 8, задача II).
Находится решение x(t ) исходной системы (3.1) с начальным
условием (3.3) и управлением u (t ) на 0; TT . На участках постоянства управления система интегрируется, как это делалось в лабораторной работе № 2. Строится график полученной линии x(t ) . Это –
траектория, исходящая из начальной точки (3.3), соответствующая
начальному сопряженному вектору (3.14) (п. 8, задача III).
Исходные данные задаются в разделе описаний констант
(см. далее текст программы, строки 4–8):
t00 – нижняя граница промежутка интегрирования;
tt – верхняя граница промежутка интегрирования;
dt1 – крупный шаг аргумента t при интегрировании систем
дифференциальных уравнений;
dt2 – мелкий шаг аргумента t при интегрировании систем дифференциальных уравнений (используется, если значение аргумента
близко к моменту переключения управления), а также шаг аргумента t при построении графика;
xmin, xmax – нижняя и верхняя границы изменения переменной x1 ;
ymin, ymax – нижняя и верхняя границы изменения переменной x 2 ;
n – количество траекторий;
alpha1 – значение угла наклона к оси абсцисс начального сопряженного вектора  0 в случае, когда нужно построить только одну траекторию, т. е. n  1 ;
grpath – константа строкового типа, путь к библиотеке графических подпрограмм, в приведенном тексте
grpath = ’d:\ turbo’;
40
41
a – двумерный массив элементов матрицы А (матрицы коэффициентов при фазовых координатах x1 и x 2 в правой части системы);
b – вектор-столбец коэффициентов при управляющем параметре u в правой части системы;
x0 – вектор-столбец начальных условий.
Результаты решения
I. Результаты, связанные с вычислениями
1. Исходные данные: матрица А (массив а), вектор b (массив b)
и начальный вектор х0 (массив х0) (см. строки 247–253).
2. Число rr и собственные значения матрицы А (см. строки
255, 256).
3. Для каждой траектории выводится номер траектории L (см.
строку 261); затем для каждого интервала постоянства управления
выводится время t – начало интервала; управление u в момент t ;
составляющие вектора состояния объекта x в момент t и составляющие сопряженного вектора P в момент t (см. строки 139–141;
148, 175, 182). Аналогичные значения t , u, x, P выводятся и для
t  TT .
II. Результаты, связанные с построением графика
Строятся последовательно N траекторий ( L  1, 2, ..., N . ). Каждая траектория окрашивается своим цветом (от «0» до «15»). Если
N  16 , то цвета траекторий (от «0» до «15») повторяются для траекторий с номерами L от 17 до 32 и т. д. (это предусмотрено в программе, строка 224).
Текст прогаммы lr3 для варианта 0
(t0 = 0, tt = 1, dt1 = 0.1, dt2 = 0.01, n = 3)
1)
2)
3)
4)
5)
6)
7)
8)
program lr3; {Траектория: dx/dt =Ax+bu, x(t00)=x00, |u|<=1}
uses CRT, Graph;
type matr = array [1..2,1..2] of real; vect = array [1..2] of real;
const t00 = 0.0; tt =1.0; ;{tt=0.931} dt1 = 0.1;dt2 =0.01;
xmin = – 3; xmax = 3; ymin = – 3; ymax = 3;
n=3; alpha1=0;{n=1; alpha1=4.48871} grPath=’d : \ turbo’ ;
e : matr = ((1,0),(0,1)); a: matr = ((–2,1),(-3, 0));
b : vect = (1, –1); x00 : vect = (2,2);
42
*
*
*
*
9)
10)
11)
12)
13)
14)
15)
16)
17)
18)
19)
20)
21)
22)
23)
24)
25)
26)
27)
28)
29)
30)
31)
32)
33)
34)
35)
36)
37)
38)
39)
40)
41)
42)
43)
44)
45)
46)
47)
48)
49)
50)
51)
52)
53)
var rr,i, j, l,u,u0,dv, mv, xc, yc : integer ;
aa, bb, t,t0,dt, kx, ky, alpha : real; ll, x, x0, x1, v, v1, p, p0, p1: vect;
function det (a : matr) : real;
begin
det : = a[1,1]*a[2,2] –a[1,2]*a[2,1]
end;
procedure eigenvals(a : matr ; var ll : vect; var aa, bb: real;
var rr: integer);
{Нахождение собственных значений матрицы A}
var p, q, d : real;
begin
p : = – ( a[1,1] + a[2,2] ); q : = det(a);
d : = sqr(p)/4 – q; aa : = – p/2;
if d > 1E–6 then
begin
bb : = sqrt (d);
ll[1] : = aa +bb ; ll[2] : = aa–bb;
rr : = 1
end
else if abs(d) < 1E–6 then
begin
ll[1] : = aa; ll[2] : = aa;
rr : = 2
end
else
begin
bb : = sqrt(–d);
rr : = 3
end
end; { eigenvals}
procedure cmatr(c : real; m: matr ; var cm : matr);
var i, j : integer ;
{Умножение матрицы на число}
begin
for i : = 1 to 2 do
for j : = 1 to 2 do
cm[i, j] : = c*m[i, j]
end;
procedure smatr (m1, m2 : matr ; var m : matr);
{Сложение матриц}
var i, j : integer ;
begin
for i : = 1 to 2 do
for j : = 1 to 2 do
m[i,j] : = m1[i,j] + m2[i,j]
43
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
54)
55)
56)
57)
58)
59)
60)
61)
62)
63)
64)
65)
66)
67)
68)
69)
70)
71)
72)
73)
74)
75)
76)
77)
78)
79)
80)
81)
82)
83)
84)
85)
86)
87)
88)
89)
90)
91)
92)
93)
94)
95)
96)
97)
98)
end;
procedure pmatr (m: matr ; v : vect ; var p : vect);
{Умножение квадратной матрицы на матрицу-столбец}
var i, j : integer ; s : real ;
begin
for i : = 1 to 2 do
begin
s : = 0;
for j : = 1 to 2 do
s : = s + m[i,j]*v[j] ;
p[i] : = s
end
end;
procedure cvect(c : real; v: vect ; var cv : vect);
{Умножение вектора на число}
var i : integer ;
begin
for i : = 1 to 2 do cv[i] : = c*v[i]
end;
procedure svect(v1, v2 : vect ; var v : vect);
{Сложение векторов}
var i : integer ;
begin
for i : = 1 to 2 do v[i] : = v1[i] + v2[i]
end;
procedure Coshy(a: matr ; v, x0: vect ; rr : integer ; t0,t : real ;
var x : vect);
{Решение задачи Коши; результат x – вектор
составляющих частного решения в момент t}
var ll1e, ll2e, aae, q1, q2, e1, e2 ,e3, expa, a1, a2 : matr;
c, c1, x0c, xc : vect;
begin
a1[1,1] : = a[2,2]; a1[1,2] : = –a[1,2];
a1[2,1] : = –a[2,1]; a1[2,2] : = a[1,1];
cmatr(1/det(a), a1, a2);
pmatr(a2, v, c);
svect(x0,c, x0c);
case rr of
1: begin
cmatr(–ll[1], e, ll1e); cmatr (–ll[2], e, ll2e);
smatr (a, ll1e, q1); smatr (a, ll2e, q2);
cmatr (exp(ll[1]*(t–t0))/(ll[1] –ll[2]), q2, e1);
cmatr (exp(ll[2]*(t–t0))/(ll[2] –ll[1]), q1, e2);
smatr (e1, e2, expa)
end;
44
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
99)
100)
101)
102)
103)
104)
105)
106)
107)
108)
109)
110)
111)
112)
113)
114)
115)
116)
117)
118)
119)
120)
121)
122)
123)
124)
125)
126)
127)
128)
129)
130)
131)
132)
133)
134)
135)
136)
137)
138)
139)
140)
141)
142)
143)
2: begin
cmatr (–ll[1], e, ll1e); smatr (a, ll1e, q1);
cmatr (t–t0, q1, e1); smatr (e, e1, e2);
cmatr(exp(ll[1]*(t–t0)), e2, expa)
end;
3: begin
cmatr (cos(bb*( t–t0)), e, e1);
cmatr (–aa, e, aae); smatr (a, aae, q1);
cmatr (sin(bb*( t–t0))/bb, q1, e2);
smatr (e1, e2, e3);
cmatr (exp(aa*(t–t0)), e3, expa)
end
end; {case}
pmatr (expa, x0c, xc);
cvect(–1,c, c1);
svect(xc, c1,x)
end; { Coshy }
function contr (p: vect): integer;
{Управление}
var i: integer; s: real;
begin
s: = 0;
for i : = 1 to 2 do s : = s + b[i]*p[i];
if s > 0 then contr: = 1 else contr: = – 1
end;{contr}
procedure intsyst (t: real; var x, p: vect);
{Интегрирование исходной и сопряженной систем}
var i, j: integer; a1: matr; v1: vect;
begin
for i : = 1 to 2 do v[i] : = b[i]*u;
Coshy(a, v, x0, rr, t0, t, x)
for i : = 1 to 2 do
begin
v1[i] : = 0;
for j : = 1 to 2 do a1[i, j] : = – a[j, i]
end;
ll[1] : = – ll[1]; ll[2] : = – ll[2];
aa : = – aa; bb : = – bb;
Coshy(a1, v1, p0, rr, t0, t, p);
ll[1] : = – ll[1]; ll[2] : = – ll[2];
aa : = – aa; bb : = – bb;
end;{intsyst}
procedure outresult(t: real; u: integer; x, p: vect);
var i : integer;
begin
45
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
144)
145)
146)
147)
148)
149)
150)
151)
152)
153)
154)
155)
156)
157)
158)
159)
160)
161)
162)
163)
164)
165)
166)
167)
168)
169)
170)
171)
172)
173)
174)
175)
176)
177)
178)
179)
180)
181)
182)
183)
184)
185)
186)
187)
188)
write (‘ t = ’ , t : 8 : 4, ‘ u = ’ , u : 2);
write (‘ x : ’ ); for i : = 1 to 2 do write (x[i] : 8 : 4);
write (‘ p : ’ ); for i : = 1 to 2 do write (p[i] : 8 : 4);
writeln
end;{outresult}
procedure traject;
begin
t0 : = t00; u0: = contr(p0);
for i : = 1 to 2 do x0[i] : = x00[i];
outresult(t0, u0, x0, p0);
for i : = 1 to 2 do
begin{Запоминание векторов x0, p0}
x1[i] : = x0[i];
p1[i] : = p0[i]
end;
dt : = dt1; t : = t0 + dt; u : = u0;
while t < tt do
begin {Основной блок}
intsyst (t, x, p); u : = contr (p);
if u = u0 then for i : = 1 to 2 do
begin {Запоминание векторов x, p}
x1[i] : = x[i];
p1[i] : = p[i]
end
else if abs(dt –dt1) < 1E–6 then
begin {Мелкий шаг}
t : = t – dt; dt : = dt2; u : = u0
end
else
begin {Новый интервал}
t0 : = t; u0 : = u; for i : = 1 to 2 do
begin
x0[i] : = x[i];
p0[i] : = p[i]
end;
write(‘ ’);
outresult (t0,u0, x0, p0);
dt : = dt1; t : = t0
end;
t : = t + dt
end; {Основной блок}
intsyst (tt, x, p);
write(‘ ’);
outresult (tt, u0, x, p)
end; {traject}
46
189)
190)
191)
192)
193)
194)
195)
196)
197)
198)
199)
200)
201)
202)
203)
204)
205)
206)
207)
208)
209)
210)
211)
212)
213)
214)
215)
216)
217)
218)
219)
220)
221)
222)
223)
224)
225)
226)
227)
228)
229)
230)
231)
232)
233)
{Далее следуют описания процедур и функций,
*
связанных с использованием графического режима}
*
*
function xscr (x: real): integer ;
{Переход от реальной (математической) координаты x к экранной} *
begin
*
xscr : = round( xc + x*kx)
*
end;
*
function yscr (y: real): integer ;
*
{Переход от реальной (математической) координаты y к экранной} *
begin
*
yscr : = round( yc– y*ky)
*
end;
*
procedure myMoveTo(x,y: real);
*
{Перемещение графического указателя
*
в точку с реальными координатами x, y}
*
begin
*
MoveTo(xscr (x), yscr (y))
*
end;
*
procedure myLineTo(x,y: real);
*
{Построение отрезка из точки, в которой находится
*
графический указатель, в точку с реальными координатами x, y}
*
begin
*
LineTo(xscr (x), yscr (y))
*
end;
*
procedure OutStr (x: real);
*
{Вывод на экран значения переменной x в формате x : 4: 1}
*
var st : string ;
*
begin
*
Str(x:4:1, st); OutText(st)
*
end;
*
procedure osy;
*
{Построение осей координат и вывод значений
*
xmin, xmax, ymin, ymax около осей}
*
begin
*
kx : = (GetMaxX – 70)/(xmax – xmin);
*
ky : = (GetMaxY – 70)/(ymax – ymin);
*
xc : = round(abs(xmin*kx)) + 30;
*
yc : = round(abs(ymax*ky)) + 30;
*
myMoveTo(xmin, 0); myLineTo(xmax,0);
*
OutTextXY(GetX, GetY – 3, ‘ >x1’);
*
myMoveTo(0, ymin); myLineTo(0, ymax);
*
OutTextXY(GetX – 3, GetY, ‘ ^x2’);
*
myMoveTo (xmin , – 0.1); OutStr(xmin);
*
myMoveTo (xmax – 0.2 , – 0.1); OutStr(xmax);
*
myMoveTo ( – 0,4, ymin + 0.1); OutStr(ymin);
*
47
234)
235)
236)
237)
238)
239)
240)
241)
242)
243)
244)
245)
246)
247)
248)
249)
250)
251)
252)
253)
254)
255)
256)
257)
258)
259)
260)
261)
262)
263)
264)
265)
266)
267)
268)
269)
270)
271)
272)
273)
274)
275)
276)
277)
278)
myMoveTo ( – 0,4, ymax – 0.1); OutStr(ymax)
end; {osy}
procedure graphic;
begin
t0 : = t00; u0: = contr (p0);
for i : = 1 to 2 do x0[i] : = x00[i];
SetColor ( l mod 16 + 1);
t : = t0;
myMoveTo(x0[1], x0[2]);
dt : = dt2; u : = u0;
while t < tt do
begin {Основной блок}
intsyst (t, x, p); u : = contr (p);
if u = u0 then
myLineTo(x[1], x[2])
else
begin {Новый интервал}
t0 : = t; u0 : = u; for i : = 1 to 2 do
begin
x0[i] : = x[i];
p0[i] : = p[i]
end
end;
t : = t +dt
end; {Основной блок}
end; {graphic}
begin {Основная программа}
ClrScr ;
{Вычисления}
writeln(‘ Матрица A
Вектор b
Вектор x00’);
for i : = 1 to 2 do
begin
for j : = 1 to 2 do write (a[i, j]:10:4);
write (b[i]:15:4, x00[i]:15:4);
writeln
end;
eigenvals(a, ll, aa, bb, rr);
writeln (‘ rr= ‘, rr);
if rr < 3 then
writeln(‘ ll[1] =’, ll[1]:8:4,‘ ll[2] = ’, ll[2]:8:4)
else
writeln(‘ ll[1] =’, aa:8:4, ‘ + ’, bb:8:4, ‘ i ll[2]= ’ ,
aa:8:4, ‘ – ’, bb:8:4, ‘ i ‘);
if n = 1 then alpha : = alpha1 else alpha : = 2*pi / n ;
for l : = 1 to n do
48
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
begin
279)
280)
TextColor(3);
281)
write(‘l= ’ , l : 2);
282)
TextColor(15);
283)
p0[1] : = cos (alpha *l); p0[2] : = sin (alpha *l);
284)
traject
end; {l}
285)
readln;
286)
287) {Построение графика}
288)
dv : = detect;
289)
InitGraph (dv, mv, grPath);
if GraphResult < > 0 then Halt(1);
290)
291)
SetBkColor(1);
292)
osy;
293)
for l : = 1 to n do
begin
294)
295)
p0[1] : = cos (alpha *l); p0[2] : = sin (alpha *l);
296)
graphic
end; {l}
297)
readln;
298)
299)
CloseGraph
300) end.
Результаты решения
I. Результаты, связанные с вычислениями
Матрица А
Вектор b
– 2.0000
1.0000
1.0000
– 3.0000
0.0000
–1.0000
*
*
*
*
*
*
*
*
*
*
Вектор x00
2.0000
2.0000
ll1  1.0000  1.4142i
ll 2  1.0000  1.4142i
aa  1.0000
bb  1.0000
rr  3
l  1 t  0.0000 u  1 x : 2.0000 2.0000 p : 0.5000 0.8660
t  0.4800 u  1
x : 0.7697 0.5452 p : 0.8765 0.8267
t  1.0000 u  1
x : 0.4567  1.0010 p : 3.7715  0.3278
l  2 t  0.0000 u  1 x : 2.0000 2.0000 p : 0.5000  0.8660
t  0.9000 u  1 x : 1.8913 1.3839 p : 0.8469  0.8058
t  1.0000 u  1 x : 0.0060 0.2804
p :  6.0940 2.2264
l  3 t  0.0000 u  1 x : 2.0000 2.0000 p : 1.0000
0.0000
t  1.0000 u  1
x : 0.1431  2.3128 p : 2.3225  1.8986
49
II. Результаты, связанные с построением графика (рис. 3.2)
Затем решаем задачи по программе lr3 с новыми значениями
N и TT , аналогичными значениям (*). Значения N и TT увеличиваются до тех пор, пока множество точек, образующих N траекторий, не будет содержать начало координат. При решении задачи варианта 0 получены значения N  150 и TT  1.5 . График траекторий
представлен на рис. 3.4.
Рис. 3.2
После решения задачи для числа траекторий N  3 строятся
графики управлений, соответствующих каждой из этих траекторий
(рис. 3.3):
l 1
l2
Рис. 3.4
После выполнения лабораторной работы студенты сдают отчет. Отчет должен содержать следующие разделы:
1. Постановка задачи и проверка условия общности положения.
2. Результаты решения задачи по программе lr3 (на языке Паскаль) для числа траекторий N  3 и TT  1, а также графики управлений u(t) для каждой из трех траекторий.
3. Значения N и TT , при которых множество точек фазовой
плоскости, образующих N траекторий, содержит начало координат
и график N траекторий для этого случая (найденные значения N
и TT будут использоваться в лабораторной работе № 4).
l 3
Рис. 3.3
50
51
Окончание таблицы
Варианты для выполнения лабораторных работ № 3, 4.
№
варианта
1
Система дифференциальных
уравнений
2
 x1  2 x1  u

 x 2  x1  5 x2  u
3
 x1   x1  3x2  u

 x 2   x2  u
x1 0   2
x 2 0   2
4
 x1  2 x1  x2  u

 x 2   x1  2 x2  u
x1 0   2
x 2 0   2
5
 x1   x1  x2  u

 x 2   x2  u
x1 0   2
x 2 0   2
6
 x1  2 x1  4 x2  u

 x 2  x1  2 x2  u
x1 0   1
x 2 0   2
7
x1   x1  u



x
 2  x1  2 x2  2u
8
 x1  3x1  x2  u

 x 2  2 x1  x2  u
x1 0   2
x 2 0   2
9
 x1   x1  x2  u

 x 2   x1  3x2  2u
x1 0   1
x 2 0    1
10
 x1   x1  x2  u

 x 2  2 x2  2u
11
 x1  2 x1  u

 x 2  2 x1  x2  u
12
 x1  5 x1  x2  u

 x 2  2 x2  u
13
 x1   x1  u

 x 2  3x1  x2  u
x1 0   2
x 2 0   2
14
 x1  2 x1  x2  u

 x 2  x1  2 x2  u
x1 0   2
x 2 0   2
 x1   x1  2 x2  u

 x 2  2 x2  u
52
Начальные условия
x1 0   2
x 2 0   2
x1 0   2
x 2 0   2
x1 0   2
x 2 0   2
№
варианта
15
Система дифференциальных
уравнений
16
 x1  2 x1  x2  u

 x 2  4 x1  2 x2  u
x1 0   2
x 2 0   1
17
 x1  2 x1  x2  2u

x 2   x2  u

x1 0   2
x 2 0   2
18
 x1   x1  2 x2  u

 x 2   x1  3x2  u
x1 0   2
x 2 0   2
19
 x1  3x1  x2  2u

 x 2  x1  x2  u
x1 0    1
x 2 0   1
20
 x1  2 x1  2u

 x2  x1  x2  u
x1 0   2
x 2 0   2
 x1   x1  u

 x2  x1  x2  u
Начальные условия
x1 0   2
x 2 0   2
Замечание 3.3. Во всех вариантах 1–20 предполагается, что
управляющий параметр u удовлетворяет ограничению u  1.
x1 0   2
x 2 0   2
x1 0   2
x 2 0   2
x1 0   2
x 2 0   2
53
ЛАБОРАТОРНАЯ РАБОТА № 4
Линейная задача оптимального быстродействия. Приближенное
построение программного управления, оптимального
по быстродействию для системы двух дифференциальных
уравнений 1-го порядка со скалярным управлением u
Необходимые теоретические сведения изложены в лабораторной работе № 3.
Итак, мы рассматриваем управляемый объект с уравнениями
движения
 x1  a11 x1  a12 x 2  b1u ,
(4.1)

 x 2  a 21 x1  a 22 x 2  b 2 u ,
ограничением на управление
u 1
(4.2)
и начальными условиями
(4.3)
x1 0   x10 , x2 0   x20 .
Требуется найти управление, переводящее объект в начало координат
x1 T   0, x2 T   0
за наименьшее время T .
Как следует из принципа максимума, необходимым и достаточным условием оптимального управления является условие
u (t )  sign(b11 (t )  b2  2 (t )),
(4.4)
где t  0, TT  (предполагается, что выполнено условие общности
положения), где  (t )  (1 (t ),  2 (t )) – решение сопряженной системы
  1   a111  a21 2 ,
(4.5)

 2  a12 1  a 22  2 .

Требуется найти такой начальный сопряженный вектор
(0)   0 , что соответствующая траектория проходит через начало
координат.
Пусть  – угол наклона начального сопряженного вектора Ψ(0)
к оси абсцисс (рис. 4.1)
(4.6)
 (0)  (cos , sin )
(предполагается, что  (0)  1 ).
2
(0)
1
Рис. 4.1
Каждому значению   0,2 соответствует функция u (t ) , получаемая по формуле (4.4), и две функции x1 (t ) и x2 (t ) – решение
системы (4.10, t  0; TT  (значение TT получено в лабораторной работе № 3, TT  Tопт , где Tопт – время оптимального быстродействия).
Введем в рассмотрение функцию r () – расстояние от траектории системы (4.1), соответствующей начальному сопряженному
вектору (4.6), до начала координат, т. е. наименьшее из расстояний
от точек траектории, соответствующей значению аргумента  , до
начала координат (рис. 4.2).
х(ТТ)
Сведение задачи к задаче минимизации функции одной переменной
Рис. 4.2
54
55
r ()  min x12 (t )  x 22 (t )  x12 (t  )  x 22 (t  ),
(4.7)
t  0;TT  .
Так как r ()  0 и равенство нулю имеет место при оптимальной траектории, имеем задачу минимизации функции (4.7) на промежутке 0    2 .
Содержание программы lr4 на языке Паскаль
При выполнении лабораторной работы № 3 получены значения N и TT , которые являются исходными данными для программы lr4. Нахождение минимального значения функции (4.7) при
  0,2 производится в два этапа.
Этап 1
1. Находим крупный шаг изменения параметра  при

  0,2 :
2   0 2
 

.
N
N
2. Находим минимальное значение функции (4.7) с крупным
шагом  , т. е. сравниваются значения r() при следующих значениях аргумента  :   0 ;    ;   2 ; …;   N  2 .
Обозначим через  значение аргумента , при котором функция
(4.7) имеет минимальное значение на 0,2 с крупным шагом:
r (  )  min r ( ) (t  0; TT ) ,
  0,2.
Этап 2
Найденное значение  уточняется.
1. Находим мелкий шаг изменения параметра  при
  [  ;   ] . Этот промежуток также разбивается на N равных частей, тогда мелкий шаг будет следующим:
(  )  (  ) 2
.
 

N
N
2. Находим минимальное значение функции (4.7) с мелким
шагом  , т. е. сравниваем значения r () при следующих значениях аргумента  :
    ;   (  )  ;   (  )  2;
~ значение параметра  , при котором функОбозначим через 
ция (4.7) имеет
  [  ;   ]
(t  0; TT ).
~
После нахождения  считаем, что решение задачи закончено.
Найдены величины: значение аргумента  , соответствующее оптимальной траектории:
~;
 опт  
время оптимального быстродействия Tопт – момент времени, в ко~) .
торый r ()  r (
Числа  опт и Tопт – результаты выполнения программы lr4.
Исходные данные к программе lr4
t 00 – нижняя граница промежутка интегрирования системы
(обычно t 00  0 );
tt – верхняя граница промежутка интегрирования системы TT ;
n – число траекторий N ;
dt1 – крупный шаг аргумента t ;
dt 2 – мелкий шаг аргумента t , а также шаг аргумента при построении графика;
x min, x max – нижняя и верхняя границы изменения переменной x1 ;
y min, y max – нижняя и верхняя границы изменения переменной x2 ;
a – матрица А коэффициентов при фазовых координатах
в правых частях уравнений системы;
b – вектор b коэффициентов при управляемом параметре u
в правых частях уравнений системы;
x 00 – вектор x 0 начальных условий.
Замечание 4.1. Шаги dt1, dt 2 берутся в 10 раз меньше, чем
в программе lr3, так как выбор оптимальной траектории требует повышенной точности.
Замечание 4.2. Текст программы lr4 будет приведен далее, см.
«выполнение лабораторной работы № 4».
  (  )  3; ... ,   (  )  N     .
56
57
Результаты решения
Выполнение лабораторной работы № 4
I. Результаты, связанные с вычислениями
1. Исходные данные: матрица А, вектор b, вектор х0 (см. строки 275–281).
2. Число rr и собственные значения матрицы А (см. строки
283–288).
3. а) для каждой траектории выводятся: номер траектории (от 0
до N ); значение параметра  (от 0 до 2π) с крупным шагом
2
; значение расстояния от траектории до начала координат
 
N
(см. строку 199);
б) голубым цветом выводятся: номер траектории, для которой
расстояние до начала координат наименьшее; значение этого расстояния; значение параметра  , соответствующего этой траектории
(   , см. строку 305);
в) для каждой траектории со значениями параметра
  [  ;   ] выводятся: номер траектории; значение пара2
метра  с мелким шагом   ; значение расстояния от траектоN
рии до начала координат (строка 199);
г) голубым цветом выводятся номер траектории с наименьшим
расстоянием, значение этого расстояния и значение параметра 
(строка 305);
д) зеленым цветом выводятся: значение tопт (время оптимального быстродействия Tопт ); значение alpha_opt (  опт , строки 315–316).
Это – результаты решения по программе lr4, их нужно записать в
тетрадь и использовать для продолжения выполнения лабораторной
работы № 4.
II. Результаты, связанные с построением графика
Выводятся графики траекторий, полученных на этапе 2, с номерами от «0» до номера оптимальной траектории (включительно).
58
Вариант 0
1. Постановка задачи
Пусть движение управляемого объекта описывается системой
дифференциальных уравнений
 x1   2 x1  x 2  u ,

 x 2   3 x1  u .
На управляющий параметр u наложено ограничение
u  1.
Начальное состояние объекта
x1 (0)  2, x2 (0)  2 .
Требуется найти управление, переводящее объект в начало координат:
x1 (T )  2, x2 (T )  2
за наименьшее время T . (Требуется приближенно построить управление, оптимальное по быстродействию.)
2. Решение задачи по программе lr4 на языке Паскаль
Замечание 4.3. Программа lr4 составлена на основе программы lr3. Для удобства строки, совпадающие со строками программы
lr3, помечены знаком «*».
Текст программы lr4 для варианта 0
(t00=0, tt=1.5, dt1=0.01,dt2=0.001, n=150)
1)
program lr4; {Оптимизация: dx/dt =Ax+bu, x(t00)=x00, x(T)=0, |u|<=1}
2)
uses CRT, Graph;
*
3) type matr = array [1..2,1..2] of real; vect = array [1..2] of real;
*
4) label 1, 2, 3;
5) const t00 = 0.0; tt =1.5; dt1 = 0.01; dt2 =0.001;
6) xmin = – 3; xmax = 3; ymin = – 3; ymax = 3;
*
7) n=150; grPath = ‘d : \ turbo’;
8) e : matr = ((1,0),(0,1)); a: matr = ((–2, 1),( –3, 0));
*
9) b: vect = (1, –1); x00 : vect = (2, 2);
*
10) var rr, i, j, l,u,u0,dv, mv, xc, yc, n_min, l0 : integer ;
11) aa, bb, t,t0,dt, kx, ky, d_alpha, alpha_1,
12) alpha_2, r, r0, rmin, t_rmin: real;
13) ll, x, x0,x1,v,v1,p, p0, p1 : vect;
14) alpha, r_alpha, t_alpha : array [0..n] of real;
*
15) function det(a : matr) : real;
*
16) begin
17)
det : = a[1,1]*a[2,2] – a[1,2]*a[2,1]
*
59
18)
19)
20)
21)
22)
23)
24)
25)
26)
27)
28)
29)
30)
31)
32)
33)
34)
35)
36)
37)
38)
39)
40)
41)
42)
43)
44)
45)
46)
47)
48)
49)
50)
51)
52)
53)
54)
55)
56)
57)
58)
59)
60)
61)
62)
end;
procedure eigenvals(a : matr ; var ll : vect; var aa, bb: real;
var rr : integer);
{Нахождение собственных значений матрицы А}
var p, q, d : real;
begin
p : = – ( a[1,1] + a[2,2] ); q : = det(a);
d : = sqr(p)/4 – q; aa : = – p/2;
if D > 1E–6 then
begin
bb : = sqrt(d);
ll[1] : = aa + bb ; ll[2] : = aa–bb;
rr : = 1
end
else if abs(d) < 1E–6 then
begin
ll[1] : = aa; ll[2] : = aa;
rr : = 2
end
else
begin
bb : = sqrt(– d);
rr : = 3
end
end; { eigenvals}
procedure cmatr(c : real; m: matr ; var cm : matr);
var i, j : integer ;
{Умножение матрицы на число}
begin
for i : = 1 to 2 do
for j : = 1 to 2 do
cm[i, j] : = c*m[i, j]
end;
procedure smatr(m1, m2 : matr ; var m : matr);
{Сложение матриц}
var i, j : integer ;
begin
for i : = 1 to 2 do
for j : = 1 to 2 do
m[i, j] : = m1[i, j] + m2[i, j]
end;
procedure pmatr (m : matr ; v : vect ; var p : vect);
{Умножение квадратной матрицы на матрицу-столбец}
var i, j : integer ; s : real ;
begin
60
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
for i : = 1 to 2 do
63)
begin
64)
65)
s : = 0;
66)
for j : = 1 to 2 do
67)
s : = s + m[i, j] * v[j] ;
68)
p[i] : = s
end
69)
70) end;
71) procedure cvect (c : real; v: vect ; var cv : vect);
72) {Умножение вектора на число}
73) var i : integer ;
74) begin
75)
for i : = 1 to 2 do cv[i] : = c*v[i]
76) end
77) procedure svect (v1, v2 : vect ; var v : vect);
78) {Сложение векторов}
79) var i : integer ;
80) begin
for i : = 1 to 2 do v[i] : = v1[i] + v2[i]
81)
82) end;
83) procedure Coshy (a: matr ; v, x0: vect ; rr : integer ; t0, t : real ;
84) var x : vect);
85) {Решение задачи Коши; результат x – вектор
86) составляющих частного решения в момент t}
87) var ll1e, ll2e, aae, q1, q2, e1, e2 ,e3, expa, a1, a2 : matr;
88) c, c1, x0c, xc : vect;
89) begin
90)
a1[1,1] : = a[2,2]; a1[1,2] : = –a[1,2];
91)
a1[2,1] : = –a[2,1]; a1[2,2] : = a[1,1];
92)
cmatr (1/det(a), a1, a2);
93)
pmatr (a2, v, c);
94)
svect (x0, c, x0c);
case rr of
95)
96)
1: begin
97)
cmatr (–ll[1], e, ll1e); cmatr (–ll[2], e, ll2e);
98)
smart (a, ll1e, q1); smart (a, ll2e, q2);
99)
cmatr (exp (ll[1]*(t–t0))/(ll[1] – ll[2]), q2, e1);
100)
cmatr (exp (ll[2]*(t–t0))/(ll[2] – ll[1]), q1, e2);
101)
smart (e1, e2, expa)
end;
102)
103)
2: begin
104)
cmatr (–ll[1],e, ll1e); smart (a, ll1e, q1);
105)
cmatr (t–t0, q1, e1); smart (e, e1, e2);
106)
cmatr (exp (ll[1]*( t–t0)), e2, expa)
end;
107)
61
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
108)
109)
110)
111)
112)
113)
114)
115)
116)
117)
118)
119)
120)
121)
122)
123)
124)
125)
126)
127)
128)
129)
130)
131)
132)
133)
134)
135)
136)
137)
138)
139)
140)
141)
142)
143)
144)
145)
146)
147)
148)
149)
150)
151)
152)
153)
154)
3: begin
cmatr (cos(bb*( t –t0)), e, e1);
cmatr (– aa, e, aae); smart (a, aae, q1);
cmatr (sin (bb* ( t–t0))/bb, q1, e2);
smart (e1, e2, e3);
cmatr (exp(aa*( t–t0)), e3, expa)
end
end; {case}
pmatr (expa, x0c, xc);
cvect (–1, c, c1);
svect(xc, c1, x)
end; { Coshy }
function contr (p: vect): integer;
{Управление}
var i: integer; s: real;
begin
s: = 0;
for i : = 1 to 2 do s : = s + b[i]*p[i];
if s > 0 then contr : = 1 else contr : = – 1
end;{contr}
procedure intsyst (t: real; var x, p: vect);
{Интегрирование исходной и сопряженной систем}
var i, j: integer; a1: matr; v1: vect;
begin
for i : = 1 to 2 do v[i] : = b[i]*u;
Coshy(a, v, x0, rr, t0, t, x);
for i : = 1 to 2 do
begin
v1[i] : = 0;
for j : = 1 to 2 do a1[i, j] : = – a[j, i]
end;
ll[1] : = – ll[1]; ll[2] : = – ll[2];
aa : = – aa; bb : = – bb;
Coshy(a1, v1, p0, rr, t0, t, p);
ll[1] : = – ll[1]; ll[2] : = – ll[2];
aa : = – aa; bb : = – bb;
end; {intsyst}
function modul (x: vect): real;
begin
modul : = sqrt( x[1]*x[1] + x[2]*x[2] )
end;
procedure traject;
begin
t0 : = t00; u0: = contr (p0);
for i : = 1 to 2 do x0[i] : = x00[i];
rmin : = modul (x0);
for i : = 1 to 2 do
62
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
begin{Запоминание векторов x0, p0}
155)
156)
x1[i] : = x0[i];
157)
p1[i] : = p0[i]
end;
158)
159)
dt : = dt1; t : = t0 + dt; u : = u0;
while t < tt do
160)
161)
begin {Основной блок}
162)
intsyst (t, x, p); u : = contr (p);
if u = u0 then
163)
begin
164)
165)
r : = modul (x);
if r < rmin then
166)
begin
167)
168)
rmin : = r; t_rmin : = t
end;
169)
for i : = 1 to 2 do
170)
begin {Запоминание векторов x, p}
171)
172)
x1[i] : = x[i];
173)
p1[i] : = p[i]
end
174)
end
175)
else if abs(dt –dt1) < 1E–6 then
176)
177)
begin {Мелкий шаг}
178)
t : = t –dt; dt : = dt2; u : = u0
end
179)
else
180)
begin {Новый интервал}
181)
182)
t0 : = t; u0 : = u; for i : = 1 to 2 do
183)
begin
184)
x0[i] : = x[i];
185)
p0[i] : = p[i]
end;
186)
187)
dt : = dt1; t : = t0
end;
188)
189)
t : = t + dt
end; {Основной блок}
190)
191)
intsyst(tt, x, p);
192)
r : = modul(x);
if r < rmin then
193)
begin
194)
195)
rmin : = r; t_rmin : = t
196)
end;
197)
r_alpha[l] : = rmin;
198)
t_alpha[l] : = t_rmin;
writeln(l : 2, alpha[l] : 8 : 3, r_alpha [l] : 8 : 3)
199)
200) end; {traject}
201) {Далее следуют описания процедур и функций,
63
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
202) связанных с использованием графического режима}
*
*
203) function xscr (x: real): integer ;
204) {Переход от реальной (математической) координаты x к экранной} *
*
205) begin
206)
xscr : = round( xc + x*kx)
*
*
207) end;
*
208) function yscr (y: real): integer ;
209) {Переход от реальной (математической) координаты y к экранной} *
*
210) begin
211)
yscr : = round( yc– y*ky)
*
*
212) end;
*
213) procedure myMoveTo(x,y: real);
214) {Перемещение графического указателя
*
215) в точку с реальными координатами x, y}
*
*
216 ) begin
217)
MoveTo(xscr (x), yscr (y))
*
*
218) end;
*
219 ) procedure myLineTo(x,y: real);
220) {Построение отрезка из точки, в которой находится
*
221) графический указатель, в точку с реальными координатами x, y}
*
*
222) begin
223)
LineTo(xscr (x), yscr (y))
*
*
224) end;
*
225) procedure OutStr (x: real);
226) {Вывод на экран значения переменной x в формате x : 4: 1}
*
*
227) var st : string ;
*
228) begin
229)
Str(x:4:1, st); OutText(st)
*
*
230) end;
*
231) procedure osy;
232) {Построение осей координат и вывод значений
*
233) xmin, xmax, ymin, ymax около осей}
*
*
234) begin
235)
kx : = (GetMaxX – 70)/(xmax – xmin);
*
236)
ky : = (GetMaxY – 70)/(ymax – ymin);
*
237)
xc : = round(abs(xmin*kx)) + 30;
*
238)
yc : = round(abs(ymax*ky)) + 30;
*
239)
myMoveTo(xmin, 0); myLineTo(xmax,0);
*
240)
OutTextXY(GetX, GetY – 3, ‘ >x1’);
*
241)
myMoveTo(0, ymin); myLineTo(0, ymax);
*
242)
OutTextXY(GetX – 3, GetY, ‘ ^x2’);
*
243)
myMoveTo (xmin , – 0.1); OutStr(xmin);
*
244)
myMoveTo (xmax – 0.2 , – 0.1); OutStr(xmax);
*
245)
myMoveTo ( – 0,4, ymin + 0.1); OutStr(ymin);
*
246)
myMoveTo ( – 0,4, ymax – 0.1); OutStr(ymax)
*
*
247) end; {osy}
*
248) procedure graphic;
64
249)
250)
251)
252)
253)
254)
255)
256)
257)
258)
259)
260)
261)
262)
263)
264)
265)
266)
267)
268)
269)
270)
271)
272)
273)
274)
275)
276)
277)
278)
279)
280)
281)
282)
283)
284)
285)
286)
287)
288)
289)
290)
291)
292)
293)
294)
295)
begin
t0 : = t00; u0: = contr (p0);
for i : = 1 to 2 do x0[i] : = x00[i];
SetColor( l mod 16 + 1);
t : = t0;
myMoveTo(x0[1], x0[2]);
dt : = dt2; u : = u0;
while t < tt do
begin {Основной блок}
intsyst (t, x, p); u : = contr (p);
if u = u0 then
myLineTo(x[1], x[2])
else
begin {Новый интервал}
t0 : = t; u0 : = u; for i : = 1 to 2 do
begin
x0[i] : = x[i];
p0[i] : = p[i]
end
end;
t : = t +dt
end; {Основной блок}
end; {graphic}
begin {Основная программа}
ClrScr ;
{Вычисления}
writeln(‘ Матрица A
Вектор b
Вектор x00’);
for i : = 1 to 2 do
begin
for j : = 1 to 2 do write (a[i, j]:10:4);
write (b[i]:15:4, x00[i]:15:4);
writeln
end;
eigenvals(a, ll, aa, bb, rr);
writeln (‘ rr= ‘, rr);
if rr < 3 then
writeln(‘ ll[1] =’, ll[1]:8:4,‘ ll[2] = ’, ll[2]:8:4)
else
writeln(‘ ll[1] =’, aa:8:4, ‘ + ’, bb:8:4, ‘ i ll[2]= ’ ,
aa:8:4, ‘ – ’, bb:8:4, ‘ i ‘);
d_alpha : = 2 * pi / n ; alpha_1 : = 0; n_min : = 1;
1: for l : = 0 to n do
begin {l}
alpha [l] : = alpha_1 + l * d_alpha;
p0 [1] : = cos ( alpha [l] ); p0 [2] : = sin (alpha [l] );
traject
end; {l}
65
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
296)
r0 : = r_alpha [0];
for l : = 1 to n do
297)
begin
298)
if r_alpha [l] < r0 then
299)
begin
300)
301)
l0 : = l; r0 : = r_alpha [l]
end
302)
end;
303)
304) TextColor (3);
305) writeln (l0 : 2, r0 : 10 : 5, alpha [l0] : 10 : 5);
306) TextColor (15);
307) readln;
308) if n_min = 2 then goto 2;
309)
alpha_1 : = alpha[l0] – d_alpha;
310)
alpha _2 : = alpha [l0] + d_ alpha;
311)
d_ alpha : = (alpha _2 – alpha _1) / n;
312)
n_min : = 2;
goto 1;
313)
314) 2 : TextColor (Green);
writeln (‘ t opt = ’, t_ alpha [l0] : 10 : 5,
315)
316)
‘, alpha opt = ’, alpha [l0] : 10 : 5);
readln;
317)
318) {Построение графика}
319)
dv : = detect;
320)
InitGraph(dv, mv, grPath);
if GraphResult < > 0 then Halt(1);
321)
322)
SetBkColor(1);
323)
osy;
for l : = 0 to n do
324)
begin {l}
325)
326)
alpha [l] : = alpha _1 +l* d_ alpha;
327)
p0[1] : = cos(alpha [l] ); p0[2] : = sin(alpha [l] );
328)
graphic;
if l = l0 then goto 3
329)
end; {l}
330)
331) 3 : readln;
332)
CloseGraph
333) end.
Исходные данные (строки 5–9):
t 00  0
tt  1 .5
n  150
получены при выполнении
лабораторной работы № 3
66
dt1  0.01 уменьшены в 10 раз по сравнению
dt 2  0.001 с программой lr3
x min  3, x max  3 , y min  3, y max  3 ;
a : matr  ((2,1), (3,0))
b : vect  (1,1); x00 : vect  (2,2)
Результаты решения
Матрица А
Вектор b
– 2.0000
1.0000
1.0000
– 3.0000
0.0000
–1.0000
ll 1  1.0000  1.4142i
ll 2  1.0000  1.4142i
Вектор x00
2.0000
2.0000
rr  3
2
,   [0;2π] (этап 1)
N
получена наилучшая траектория с номером l 0  107 , расстоянием
до начала координат r 0  0.00604 , значением параметра  :
(l 0)  4.48201   (голубой цвет).
2
,
 
После
решения
с
мелким
шагом
N
  [  ;   ] (этап 2) получена наилучшая траектория с номером l 0  87 , расстоянием до начала координат r 0  0.00123
и значением параметра  :
(l 0)  4.48201   (голубой цвет).
Далее получается:
tопт  0.93100;
После решения с крупным шагом  
*
*
*
*
*
*
*
*
*
alpha _ opt  4.
48871 (зеленый цвет).
Эти числа следует записать в тетрадь.
После этого строятся графики траекторий, полученных на этапе 2, последняя траектория с номером l 0  87 – оптимальная.
3. Построение оптимальной траектории с помощью программы lr3 на языке Паскаль
Текст программы lr3 совпадает с текстом, приведенном в предыдущем разделе, за исключением следующих изменений:
67
tt  0.93100
n 1
alpha1  4.48871
dt1  0.01
dt 2  0.01
(строка 4),
(строка 6),
(строка 6),
(строка 4),
(строка 4).
Результаты решения
а) Матрица A
Вектор b
– 2.0000
1.0000
1.0000
– 3.0000
0.0000
–1.0000
ll 1  1.0000  1.4142i
ll 2  1.0000  1.4142i
l  1 : t  0.0000 u  1
rr  3
x : 2.00 0;
2.0000
p :  0.2200;  0.9800
t  0.1710 u  1
x : 1.7663;
0.8580
p :  0.8878;  0.8828
t  0.9310 u  1
x : 0.0003;  0.0012
p :  5.6059; 1.4562
б) строится график оптимальной траектории (рис. 4.3):
Используя результаты решения по программе lr3, строим график оптимального управления (вручную или с помощью компьютера) (рис. 4.4).
Вектор x00
2.0000
2.0000
Рис. 4.4
4. Оценим точность попадания траектории в начало координат за время t  Tопт .
Нужно найти число , такое, что
x1 (Tопт )   x1 (0) ,
x2 (Tопт )   x2 (0) .
Решение:
Рис. 4.3
 x1 (Tопт ) x 2 (Tопт ) 
  max 0.0003 ; 0.0012  
  max
;
 x ( 0)
x 2 (0) 
2 
 2
 1
 max(0.00015;0.0006)  0.0006.
После выполнения лабораторной работы № 4 студенты сдают
отчет, который включает следующие пункты:
1. Постановка задачи.
2. Результаты решения по программе lr4: числа Tопт ,  опт .
3. Результаты решения по программе lr3: траектории, моменты
переключения оптимального управления, значения фазовых координат и управляющего параметра на концах промежутков постоянства управления.
4. График оптимального управления.
5. Оценка точности попадания траектории в начало координат
за время Tопт .
68
69
Варианты заданий для выполнения лабораторной работы
№ 4 совпадают с вариантами заданий для выполнения лабораторной
работы № 3 и приводятся в конце предыдущего раздела («Лабораторная работа № 3).
70
Рекомендуемая литература
1. Понтрягин Л. С. Математическая теория оптимальных процессов /
Л. С. Понтрягин, В. Г. Болтянский, Р. В. Гамкрелидзе, Е. Ф. Мищенко. – М.:
Наука, 1969. – 384 с.
2. Болтянский В. Г. Математические методы оптимального управления
/ В. Г. Болтянский. – М.: Наука, 1969. – 408 с.
3. Баринов Н. Г. Оптимизация процессов и систем управления в судоавтоматике / Н. Г. Баринов. – Л.: Судостроение, 1976. – 256 с.
4. Эрроусмит Д. Обыкновенные дифференциальные уравнения. Качественная теория с приложениями / Д. Эрроусмит, К. Плейс. – М.: Мир, 1986. –
240 с.
5. Нарбут Л. К. Системы линейных дифференциальных уравнений
и оптимальное управление. Ч. 1. Решение задачи Коши для систем линейных
дифференциальных уравнений первого порядка с постоянными коэффициентами / Л. К. Нарбут; СПбГАСУ. – СПб., 2009. – 180 с.
6. Марченко А. И. Программирование в среде Turbo Pascal 7.0 / A. И. Марченко, Л. А. Марченко. – М.: Бином Универсал; Киев: Юниор, 1997. – 494 с.
7. Культин Н. Б. Программирование в Turbo Pascal 7.0 и Delphi /
Н. Б. Культин. – СПб.: БХВ – Петербург, 2007. – 400 с.
71
Оглавление
Предисловие……………………………………………………………….
Введение……………………………………………………………………
Лабораторная работа № 1
Решение задачи Коши для системы двух линейных однородных
дифференциальных уравнений 1-го порядка с постоянными
коэффициентами……………………………………………………..
Выполнение лабораторной работы № 1……………………………
Лабораторная работа № 2
Решение задачи Коши для простейшей системы двух линейных
неоднородных дифференциальных уравнений 1-го порядка
с постоянными коэффициентами……………………………………
Выполнение лабораторной работы № 2……………………………
Лабораторная работа № 3
Линейная задача оптимального быстродействия. Построение
траекторий для системы двух дифференциальных уравнений
1-го порядка со скалярным управлением u…………………………
Выполнение лабораторной работы № 3……………………………
Лабораторная работа № 4
Линейная задача оптимального быстродействия. Приближенное
построение программного управления, оптимального
по быстродействию для системы двух дифференциальных
уравнений 1-го порядка со скалярным управлением u…………
Выполнение лабораторной работы № 4…………………………
Рекомендуемая литература………………………………………
Учебное издание
Нарбут Людмила Константиновна
Букунова Ольга Викторовна
ЛАБОРАТОРНЫЕ РАБОТЫ ПО ТЕОРИИ УПРАВЛЕНИЯ
Учебное пособие
Редактор О. Д. Камнева
Корректор К. И. Бойкова
Компьютерная верстка И. А. Яблоковой
Подписано к печати 03.12.13. Формат 6084 1/16. Бум. офсетная.
Усл. печ. л. 4,2. Тираж 100 экз. Заказ. 188. «С» 77.
Санкт-Петербургский государственный архитектурно-строительный университет.
190005, Санкт-Петербург, 2-я Красноармейская ул., д. 4.
Отпечатано на ризографе. 190005, Санкт-Петербург, 2-я Красноармейская ул., д. 5.
72
3
4
5
7
20
21
34
38
54
59
71
Документ
Категория
Без категории
Просмотров
1
Размер файла
822 Кб
Теги
labor, rabota, narbut
1/--страниц
Пожаловаться на содержимое документа