Воронов А.А. Теория автоматического управления. Часть первая

Подождите немного. Документ загружается.

ние двух членов ряда Тейлора. Наклон решений х (t) в началь-

ной точке определяется по формуле

**f(x

Приближение х

кх (У находится с помощью двух первых

членов ряда:

ik) » = + hf (x

Полагая t

~ t

+ К находим х (*

) ^ х

= х

+ hf (x

)

и т. д. Этот процесс можно продолжить по формуле

+ У»

= 1, 2, ....

Погрешность метода имеет порядок О (А

), так как члены ря-

да Тейлора, содержащие h во второй и более высоких степенях,

отбрасываются.

Точность метода можно увеличить на порядок, если исполь-

зовать среднее значение производной в начале и конце интер-

вала шага. Геометрически это означает, что наклон касательной

на середине шага меняется. Такой усовершенствованный метод

Эйлера x

h+1

= x

+ hf (x

, t

h) согласуется с разложением

в ряд Тейлора вплоть до членов степени

Порядок ошибки со-

ставляет О (/t

Расчетная формула имеет вид

+ 0 МКг + Ы*

где

Ki*=hf(x

); K

~hf(x

+ Ki, t

+Kh £-0, 1, ..., k — l.

Таким образом, простой и усовершенствованный методы

Эйлера можно рассматривать как приближение, использую-

щее два и три члена ряда Тейлора соответственно.

Если использовать большее число членов ряда Тейлора

x{t + h)^x(t) + hx(t) + -~- x(t) + ...

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

Использование первых пяти членов ряда Тейлора дает

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

ности

4-4" Wi + + +

*> •

где К! ~f(x

, t

); Кг = f(t

+ Л/2,

+ К,/2);

=f(t

+ Л/2, + /С

/2); Кь = f (t

+ h, x

+ K

В методе Рунге—Кутта четвертого порядка не требуется

вычислений производных, функция / вычисляется четыре ра-

за для продвижения на один шаг вперед. Так как формулы опи-

сывают метод четвертого порядка точности, то порядок по-

грешности метода составляет О (Л

). Для выбора шага может

л Г*»-**!

быть использована Оценка к, , которая не должна пре-

вышать нескольких сотых. Имеются модификации формул

Рунге—Кутта. Одна из них, предложенная Мерсоном, поз-

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

точность. Используемые формулы Рунге—Кутта—Мерсона

имеют вид

**+!

^ *н + 0,5 (Кг + 4/(

+ /С

где

K^-jhf (x

, t

); K*~hf ( t

+ ~h

+ Кг);

Ks —-J hf [t

+ x

+ ± Кг + ~ K*);

Погрешность по этому методу оценивается по формуле

е ==

/с

_|./С

+ 4К*-К»,

где е — заданная точность.

Если правая часть превышает заданную погрешность бо-

лее чем в пять раз, шаг h уменьшается в два раза, если правая

часть меньше чем 5/32 заданной погрешности, то шаг удваи-

вается. Автоматическое изменение шага, по утверждению Мер-

сона, на 20% ускоряет процесс по сравнению со стандартной

процедурой Рунге—Кутта с постоянным шагом. Экономия вре-

мени достигается за счет того, что при стандартном методе

шаг выбирается слишком малым и время счета увеличивается.

Можно построить формулы Рунге—Кутта высших степеней,

при этом основная часть расчетов приходится на счет правой

части уравнения. Формула степени точности р требует р-крат-

ного вычисления правой части. Это может привести к значи-

тельному увеличению затрат машинного времени.

другой сто-

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

большего шага h. Методы Рунге—Кутта легко программиру-

ются.

Величину шага h можно изменить на любом этапе вычисле-

ний. Метод является «самостартующим», первые точки реше-

ния вычисляются естественно, как и все остальные.

Перечисленные методы могут быть как явными, так и не-

явными. Явными методы называются по той причине, что ис-

комое значение на (k + 1)-м шаге выражается явно через

значения х

и F

, полученные на предыдущих шагах. Напри-

мер, переходные процессы для явного метода Эйлера рассчиты-

вают по формуле

h+l

+ hF(t

Неявные методы — это такие, в которых искомое значение

k+i определяется неявно, т. е. через значения х

к+1

и F

на том же шаге. Для неявного метода формула Эйлера имеет

вид

—

k + hF (х

к+ъ

fft+i).

Важно, что в этом методе можно брать любой шаг, меняя лишь

точность построения процессов. Метод устойчив при любых

значениях h (h > 0).

Вычисление x

в неявных методах сложнее, чем в явных,

так как приходится решать систему алгебраических уравне-

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

интегрировании определенных классов систем общая трудоем-

кость может быть меньше.

Ко второму классу относятся многошаговые (многоступенча-

тые) методы. Их отличительная черта — использование ин-

формации при вычислении следующей точки t

) не

только в точке (x

)

но и в предыдущих точках. Многоша-

говые методы послужили базой для создания методов с прог-

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

сказывается значение х

(прогноз), а затем оно каким-либо

способом исправляется (коррекция). Для корректировки мож-

но использовать ту же самую формулу. Итерационный процесс

повторяется до тех пор, пока не достигается заданная точность.

Однако, как правило, используются две формулы, называе-

мые соответственно формулами прогноза и коррекции. Так

как многошаговые методы не обладают «стартовым» свойством,

. TP начинать решение надо с помощью одношаговых методов.

Чаще всего для начала решения используется метод Рунге—

Кутта. В настоящее время для интегрирования систем х =

= / (х, t) широко используются методы Адамса—Башфорта и

Хэмминга.

Общая формула прогноза для методов четвертого порядка

точности имеет вид

н + А х

+ Л

+ h (В

+ В

_, + В

х^

+ В

х'ь^з) +

(-^gp-),

где 4);

Л-

—Л— Л

; 59+ ЮЛ^ЗгЛа)^;

Л,, = Ац В

= (37—5А

+ 8А

)/24;

-Л

; В

-(-9+Л)/24;

(55

+ 9Л

+ 8Л

)/24; £—(—251

—

19Лх—8Л

)/6.

Формула прогноза типа Адамса—Башфорта может быть

получена интегрированием обратной интерполяционной фор-

мулы Ньютона 117, 18].

Прогноз по методу Адамса—Башфорта осуществляется

по формуле

Xk+i = ——

Л (55х£ —

59х^ -1 + 37*£

—

Коррекция выполняется по формуле

ьлг

+ — h {9x'

+ 194 — 5xfe_

+ x'k-2) +

L 720 J

Последние члены в обеих формулах в действительности в вы-

числениях не используются и служат для оценки ошибок дис-

кретизации (усечения).

В распространенном в настоящее время методе Хемминга

используются следующие формулы прогноза и коррекции:

п р о г н о з-

4 оо

- х„_з +-jh (24 -4- 1+24- ft

коррекция

h+l

=4-- +

зл

W+1 + 24—4- i)]

х<5)

о 40

Для получения требуемой точности формулы прогноза и

коррекции должны быть одного порядка. Особенность методов

с прогнозом и коррекцией состоит в том, что они позволяют

находить разность между прогнозируемым и скорректирован-

ным значениями и устранять ошибку. Многошаговые методы

более экономичны в смысле затрат машинного времени, так как

используют информацию о ранее вычисленных точках. Однако

при любом изменении величины шага h приходится временно

возвращаться к одношаговым методам. Методы, разработан-

ные в самое последнее время, позволяют менять порядок точ-

ности и шаг. В качестве корректирующей часто используется

неявная формула, в которую подставляются данные прогноза.

§ 6.8. Машинная реализация частотных методов

Машинная реализация известных, сложившихся в теории

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

пространить их на системы большой размерности, выполнять

исследования не только по одному, но и по нескольким пара-

метрам, упростить и ускорить процедуру получения конеч-

ных данных, осуществить сервисное представление результа-

тов (графиков, таблиц, расчетных данных) с помощью внеш-

них устройств ЭВМ..

Несмотря на наметившуюся тенденцию широкого внедре-

ния ЦВМ в область анализа и синтеза автоматических систем,

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

утратили своего значения. Наоборот, реализация их на ЦВМ

позволяет в кратчайший срок получать обширную и весьма

ценную информацию о проектируемой системе. По амплитудно-

фазовым частотным характеристикам (АФЧХ) проектировщик

может судить о таких качественных характеристиках, как за-

пасы устойчивости по амплитуде и по фазе, резонансная ча-

стота и т. д. Исследование реальных объектов с помощью

АФЧХ дает возможность решать задачи анализа функциональ-

ных, структурных и параметрических свойств объекта и отдель-

ных его частей, идентификации по экспериментально снятым

АФЧХ, синтеза регулятора систем путем подбора корректи-

рующих контуров.

Рассмотрим возможности вычисления на ЦВМ амплитудно-

фазовых частотных характеристик в логарифмическом масшта-

бе (ЛАФЧХ) и укажем на особенности и основные трудности,

возникающие при решении этой задачи на машине.

Пусть передаточная функция задана в виде

«"-H-.-.+frm

(б48>

Q(s)

Требуется вычислить U (со) = Re W (/со) и V (со) =

= Iт W (/со), а затем построить

А (со) - 20 lg |/(У

(со)

+ V

(со)

Ф (со)

= arctg

(сo)/U

(со))

+ гал, т = 0, ± 1, ±2,

где А (со) и ф (со) — амплитудная и фазовая характеристики

исследуемой системы.

Построение частотных характеристик сводится к много-

кратному вычислению передаточной функции W (s) при ком-

плексном значении аргумента s = I + /со, где /

= —1, а со

принимает значения из некоторого интервала l<o

, co

max

и последующему построению графиков модуля и аргумента.

Факт устойчивости или неустойчивости устанавливается по

подсчету точек пересечения фазовой характеристикой линии

= —п слева от частоты среза lg {А (со)} = 0. Машинная

ориентация этого метода оценки устойчивости состоит в алго-

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

А (со), ф (со). Особенность заключается в том, что функция

Ф (со) может находиться не только в первом квадранте, т. е.

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

при машинной реализации, ибо на ЦВМ легко могут быть по-

лучены только главные значения функции arctg х.

Рассмотрим один из методов построения ЛАФЧХ. Пред-

ставим передаточную функцию W (s) в виде комбинации эле-

ментарных функций

(S) = (а! + р? +с?) > 0,

которые соответствуют передаточным функциям типовых зве-

ньев системы. Если принять во внимание, запаздывание в бло-

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

(s)

У*

. (6.49)

aiS^

+ biS+Ci

Таким образом, можно представить передаточную функцию

любого элементарного звена, полагая соответствующие ко-

эффициенты а

, р

-, Yf, a

равными нулю. Программа

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

чение полной передаточной функции получается путем эле-

ментарных арифметических действий над комплексными числа-

ми.

При построении ЛАФЧХ на ось абсцисс, как обычно, на-

носятся значения lg со, меняющиеся в интервале [lg co

lg ю

тах

]. На оси ординат откладываются величины 20 lg А (со)

и ф (со). Программа вычисления величины 20 lg А (со) и

<р

(о)

состоит из постоянной части, которая содержит стандартную

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

ции (6.49), а также стандартной подпрограммы вычисления

функций e

, sin сот и cos сот и подпрограммы сложения и

умножения пар комплексных чисел. При проектировании кон-

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

ды и выходы соответствующих элементарных передаточных

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

ным правилам.

На ЦВМ осуществляется расчет 20 lg А

(со) и <р

(со),

i £ [1, А/], по стандартным программам, составленным в со-

ответствии с выражениями

Л, fo)

Jh-JLL-

\e"**\

(6.50)

Ci Y(1

—

(О

)

-f (b\ <0)2

P i

i <°

q>, (со)

- arctg /

—arctg , + arctg

cot,

(6.51)

I ^

где a I =a

; p/ -

/Y«;

a/ b\ - bjc^

Используя правила разложения полиномов на множители

+ + уи a

+ biS + c

исходную передаточную функ-

цию приводят к произведению из звеньев типа (6.49). Тогда

логарифмическая амплитудно-частотная характеристика ра-

зомкнутой передаточной функции W (s), 20 lg А (со) будет

иметь вид

20 lg А (со) 20 lg Ai (со),

а фазовая

ф(°>) = 2 ч^®)-

/ « I

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

вить программу вычислений ЛАФЧХ. Ввод массива коэффи-

циентов исходных типовых звеньев осуществляется после

их идентификации.

Остановимся на некоторых особенностях построения

ЛАФЧХ по приведенной методике. При реализации формулы

(6.51) для углов ф»я/2 значения аргументов в (6.51) стремят-

ся к оо, что исключает расчет таких точек на ЦВМ, поэтому их

следует заменить на выражения

[Wi (/со)}

. Im{r

(/w)}

cpi = arccos I—

или ф; = arcsin -—IILJLL

warn i(/o>)i -

которые не имеют указанных особенностей.

В общем случае при построении характеристики <р

(со)

для звеньев Wi

(s)

может быть скачок фазы при следующих соче-

таниях коэффициентов Ъи Pi'

bt= 0; Pi ф 0; Ъ,ф 0; р,-0; 0.

В связи с этим программа вычисления фазы должна стро-

иться с учетом возможных скачков функции Ф (со). Если

ф 0, bi Ф 0, т. е. система демпфирована, то ф (со) — не-

прерывная функция и скачка фазы не будет.

В любом из указанных «опасных» сочетаний коэффициентов

следует проверить, попадают ли значения

= V Усo

= Vcjot

в интервал изменения частоты (co

, co

fe+1

). Если со

/г

- е Ico

ft+1

J, /—1,2, то функция на этом интервале терпит разрыв

и фаза меняется скачком на величину ±я. Тогда к значению

Ф (co

) на этом интервале следует добавить ±я по формуле

фt (со) = фi (со) ± я.

Для учета скачка фазы можно также использовать следую-

щий подход. Пусть известно истинное значение фазы ф (со) в

точке со = Выполним вычисления для получения частот-

ных характеристик в точке co

fe+1

. Главное значение функции

<р

(<*>

Л+1

) в точке со = toft

определяется по формуле (6.51):

р;

со

ь[ о

ф* (со)

- arctg ~ —"

—

arctg

—

—- + arctg

сот.

1—a} о>

1—a

Истинное значение ф (<%+*) в случае скачка фазы мо-

отличатся отг л ав ног о значения ф (co

fe+1

) слагаемым, рав-

ным /Ш1, тде ^ drl- Будем считать, что co

fe+1

— h<x>

где

h —шаг в логарифмическом масштабе. Если шаг h достаточно

мало отличается от единицы, то можно считать, что в случае

скачка фазы разность ф (со

) — ф (со

й+1

) мало отличается от

лис. Тогда значение фазовой характеристики в точке co

fe+1

= ha)

вычисляется по формуле

ф<«*+,)-»(с) j

(652)

где Е (х) = т.

Если скачка фазы нет, то Е (х) округляется до нуля.

Шаг h построения частотной характеристики должен быть

выбран таким, чтобы колебания фазовой характеристики при

изменении частоты от со

до co

fe+1

не превышали значения зх.

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

20 lg А (со) необходимо вычислить значение амплитудной ха-

рактеристики при том значении со, при котором фазовая харак-

теристика обращается в —л. Следовательно, нужно найти ко-

рень фазовой характеристики и при этом значении корня вы-

числить значение амплитудной характеристики.

Значения co

mln

, со

со

, ..., о>

таХ

находятся в памяти маши-

ны в виде таблицы. Следовательно, решение задачи сводится

в основном к решению уравнения ф (со) = —я, левая часть

которого задана таблично.

При построении частотного годографа значения со& и co

fe+1

располагаются достаточно близко, поэтому для эффективного

решения уравнения ф (со) — —л может быть использован ли-

нейный интерполяционный подход. Использование линейной

интерполяции по формуле

т г*/,л \ <°/гН

COFC

+ Ф (СО,,)

позволяет быстро и достаточно точно определять значение,

при котором фазовая характеристика ф (со) ==. — л. ;

В том случае, когда передаточную функцию W (s) предста-

вить в виде звеньев Wi (s) не удается, приходится строить ал-

горитм вычисления ЛАФЧХ для общего случая представления

в виде (6.48).

В результате внедрения таких программ инженеру-про-

ектировщику остается следующая работа: а) написать выра-

жения для передаточных функций в виде комбинаций переда-

точных функций для элементарных звеньев; б) указать, от ка-

кого входа до какого выхода системы необходимо пострс&ть

частотные характеристики.

Использование ЦВМ для получения частотных характери-

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

ния автоматизации вычислений, но и принципиальные преиму-

щества по сравнению с обычными «ручными» способами. При

исследовании большого числа точек пространства параметров

целесообразно выводить на печать значения функций А (со)

и ф (со) не при всех

(k = 0, 1, 2, ...), а при тех, которые в

наибольшей степени интересуют исследователя. В связи с этим

необходимо выявить те параметры, которые анализирует ис-

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

ной характеристики. Обычно это запасы устойчивости по ам-

плитуде и по фазе, которые могут быть легко получены при

наличии этой характеристики, а также резонансный пик и ча-

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

А (со) и ф (со) для большого числа co

, а выводить на печать

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

ния запасов устойчивости по амплитуде и по фазе, величины

резонансной частоты и т. д.

Составленные подобным образом программы прошли до-

статочную практическую проверку и позволили в полной мере

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

плитудно-фазовых частотных характеристик на ЦВМ. Частот-

ные методы были реализованы на ЦВМ (в кодах машины) и

применялись В. М. Есиповым еще в начале 60-х годов при ма-

шинном синтезе систем управления.

Можно строить частотные характеристики путем прямого

вычисления вынужденной составляющей уравнения перемен-

ных состояния при подаче на вход гармонического сигнала

A sin со/. Такой способ построения АФЧХ был предложен

Р. И. Сольницевым. Преимущества его состоят в том, что пред-

ставляется возможным строить частотные характеристики не-

посредственно по уравнениям переменных состояния (ибо ча-

стотная характеристика, по существу, представляет собой ре-