Лазарев Ю.Ф. Mатематическое моделирование физических процессов и технических систем в MATLAB
Подождите немного. Документ загружается.
221
Zero/pole/gain from input 1 to output...
s
#1: --------------
s (s^2 + 100)
10
#2: --------------
s (s^2 + 100)
Zero/pole/gain from input 2 to output...
-10
#1: --------------
s (s^2 + 100)
s
#2: --------------
s (s^2 + 100)
Ввиду того, что первая (ss) модель GYROss была создана непосредственно процедурой-конструктором по
заданным числовым данным, а последующие (GYROtf и Gyrozp) - путем преобразования уже созданной
модели, будем называть модель, созданную конструктором, основной, а остальные - вспомогательными.
Отметим, что ss-модель в MatLAB можно создать и по системе дифференциальных уравнений первого
порядка, не разрешенных относительно
производных, т.е. когда система описывается совокупностью уравнений
вида:
⎪
⎩
⎪
⎨
⎧
⋅+⋅=
⋅+⋅=⋅
,uDxCy
uBxA
dt
dx
E
(6.9)
где - произвольная квадратная матрица размером (n×n), а n - порядок заданной системы дифференциальных
уравнений. Для этого следует уже использовать не конструктор
ss
, а специальную процедуру
dss, отличие
которой от предыдущей лишь в том, что она требует задания не четырех, а пяти матриц, последней из которых
должна быть матрица Е.
E
В качестве примера рассмотрим уравнения того же гироскопа в виде
⎪
⎩
⎪
⎨
⎧
=⋅−⋅
=⋅+⋅
).(
)(
2
1
tLHJ
tNHJ
αβ
βα
&
&&
&
&&
(6.10)
Вводя те же переменные (см. (3), (5)...(7)), получим систему уравнений в виде (9), где матрицы A и E будут
иметь вид:
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
−
=
000
000
1000
0100
H
H
A
. (6.11)
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
1
1
000
000
0010
0001
J
J
E
Остальные матрицы - B, C и D будут прежними (8). Введем новые матрицы при таких значениях параметров -
H=10, J
1
=2, J
2
=3:
H= 10; J1=2; J2=3;
A=zeros(4); A(1,3)=1; A(2,4)=1; A(3,4)=-H; A(4,3)=H
A =
0 0 1 0
0 0 0 1
0 0 0 -10
0 0 10 0
E=eye(4); E(3,3)=J1; E(4,4)=J2
E =
222
1 0 0 0
0 1 0 0
0 0 2 0
0 0 0 3
Теперь зададим ss-модель, пользуясь процедурой dss:
Gyross=dss(A,B,C,0,E)
a =
x1 x2 x3 x4
x1 0 0 1 0
x2 0 0 0 1
x3 0 0 0 -10
x4 0 0 10 0
b =
u1 u2
x1 0 0
x2 0 0
x3 1 0
x4 0 1
c =
x1 x2 x3 x4
y1 1 0 0 0
y2 0 1 0 0
d =
u1 u2
y1 0 0
y2 0 0
e =
x1 x2 x3 x4
x1 1 0 0 0
x2 0 1 0 0
x3 0 0 2 0
x4 0 0 0 3
Continuous-time model.Continuous-time system.
Как и ранее, создадим на этой основе вспомогательные tf- и zpk- модели:
Gyrotf=tf(Gyross)
Transfer function from input 1 to output...
0.5 s - 1.11e-016
#1: -----------------
s^3 + 16.67 s
1.667
#2: -------------
s^3 + 16.67 s
Transfer function from input 2 to output...
-1.667
#1: -------------
s^3 + 16.67 s
0.3333 s + 1.48e-016
#2: --------------------
s^3 + 16.67 s
Gyrozp=zpk(Gyross)
Zero/pole/gain from input 1 to output...
0.5 s
#1: ----------------
s (s^2 + 16.67)
1.6667
223
#2: ----------------
s (s^2 + 16.67)
Zero/pole/gain from input 2 to output...
-1.6667
#1: ----------------
s (s^2 + 16.67)
0.33333 s
#2: ----------------
s (s^2 + 16.67)
В предыдущих примерах за основу была принята ss-модель. Но в качестве основной можно выбрать и любую
из двух других моделей. Примем, например, в качестве основной модель в передаточных функциях.
Система, описываемая уравнениями (10), если принять те же, что и ранее входные и выходные величины, имеет
4 передаточных функции, которые образуют матрицу
передаточных функций размером (2×2). Каждый из
столбцов этой матрицы содержит передаточные функции, соответствующие некоторой одной входной
величине по всем выходным величинам. Определенная строка матрицы, наоборот, содержит передаточные
функции какой-то одной выходной величины по всем входам системы. В целом матрица передаточных
функций в рассматриваемом случае может быть представлена в виде:
⎥
⎦
⎤
⎢
⎣
⎡
=
)()(
)()(
)(
2221
1211
sWsW
sWsW
sW
В соответствии с уравнениями (10) значения элементов этой матрицы равны:
;)(
22
21
2
11
HsJJ
J
sW
+⋅⋅
=
;
)(
)(
22
21
12
HsJJs
H
sW
+⋅⋅⋅
−=
;
)(
)(
22
21
21
HsJJs
H
sW
+⋅⋅⋅
=
22
21
1
22
)(
HsJJ
J
sW
+⋅⋅
=
. (6.12)
Чтобы ввести эти передаточные функции и создать на их основе tf-модель, следует вначале создать два массива
ячеек
- массив ячеек размером (2×2) из векторов коэффициентов всех числителей передаточных функций;
- массив ячеек такого же размера из векторов коэффициентов знаменателей передаточных функций.
Для рассматриваемого случая это можно сделать так. Сначала создадим вектор
коэффициентов общей части
знаменателей:
[
]
2
211
,0, HJJV
zn
⋅=
,
затем - вектор дополнительного множителя в некоторых знаменателях:
[]
0,1=V
.
Теперь создадим вектор коэффициентов второго знаменателя путем свертки этих двух векторов (это
соответствует перемножению полиномов):
),(
12
VVconvV
znzn
=
.
Сформируем массив den ячеек знаменателей по схеме:
for k1=1:4
for k2=1:2
den(k1,k2)={
V
};
end
end
den(1,2)={
V
}; den(2,1)={
V
}.
zn1
zn2 zn2
224
Переходя к определению массива ячеек nom числителя, можно записать его таким образом: nom = {J
2
, -H; H,
J
1
}.
Теперь можно сформировать tf-модель, используя установленные матрицы ячеек числителей и знаменателей.
Ниже приводится пример этого:
>> Vzn1=[J1*J2, 0, H^2]
Vzn1 =
6 0 100
>> V=[1, 0]
V = 1 0
>> Vzn2=conv(Vzn1,V)
Vzn2 = 6 0 100 0
>> for k1=1:2
for k2=1:2
den(k1,k2)={Vzn1};
end
end
>> den(1,2)={Vzn2}; den(2,1)={Vzn2}
den =
[1x3 double] [1x4 double]
[1x4 double] [1x3 double]
>> nom ={J2, -H; H, J1}
nom =
[ 3] [-10]
[10] [ 2]
>> gyrotf=tf(nom,den)
Transfer function from input 1 to output...
3
#1: -----------
6 s^2 + 100
10
#2: -------------
6 s^3 + 100 s
Transfer function from input 2 to output...
-10
#1: -------------
6 s^3 + 100 s
2
#2: -----------
6 s^2 + 100
225
Предостережение. При манипуляциях или преобразованиях LTI-модели следует учитывать, что:
1) три формы представления LTI-объектов не являются эквивалентными при
численных расчетах; в частности, точность вычислений с передаточными функциями
высокого порядка часто недостаточна; старайтесь работать по преимуществу со
сбалансированными моделями пространства состояния и использовать передаточные
функции только для отображения на экране или для
интерпретации (расшифровки)
результатов;
2) преобразования в формат передаточных функций может сопровождаться потерями
точности; в результате, полюсы передаточной функции могут заметно отличаться от
полюсов заданной zpк-модели или модели пространства состояния (для проверки
наберите
help roots);
3) преобразования в пространство состояния не являются однозначно определенными
в случае одномерной системы и не гарантируют создания минимальной
конфигурации системы в случае многомерной системы; так, заданная в пространстве
состояния модель sys, при преобразовании
ss(tf (sys)) может сформировать
модель с другими матрицами пространства состояния или даже с другим числом
переменных состояния в многомерном случае; таким образом,
следует по
возможности избегать преобразований моделей из одной формы в другую и
наоборот.
Проиллюстрируем это на примере. Преобразуем созданную основную TF-модель gyrotf в SS-модель gyross:
gyross=ss(gyrotf)
a =
x1 x2 x3 x4
x1 0 -4.167 0 0
x2 4 0 0 0
x3 0 0 0 -4.167
x4 0 0 4 0
x5 0 0 0 4
x6 0 0 0 0
x7 0 0 0 0
x8 0 0 0 0
x9 0 0 0 0
x10 0 0 0 0
x5 x6 x7 x8
x1 0 0 0 0
x2 0 0 0 0
x3 0 0 0 0
x4 0 0 0 0
x5 0 0 0 0
x6 0 0 -4.167 0
x7 0 4 0 0
x8 0 0 4 0
x9 0 0 0 0
x10 0 0 0 0
x9 x10
x1 0 0
x2 0 0
x3 0 0
x4 0 0
x5 0 0
x6 0 0
x7 0 0
x8 0 0
x9 0 -4.167
x10 4 0
b =
u1 u2
226
x1 0.25 0
x2 0 0
x3 0.25 0
x4 0 0
x5 0 0
x6 0 0.25
x7 0 0
x8 0 0
x9 0 0.25
x10 0 0
c =
x1 x2 x3
y1 0 0.5 0
y2 0 0 0
x4 x5 x6
y1 0 0 0
y2 0 0.4167 0
x7 x8 x9
y1 0 -0.4167 0
y2 0 0 0
x10
y1 0
y2 0.3333
d =
u1 u2
y1 0 0
y2 0 0
Continuous-time model.
Нетрудно убедиться, что мы получили систему, совсем не похожую на ранее введенную S-модель (11), хотя обе
они описывают одну и ту же ЛСС. Новая модель отличается от предыдущей не только иными значениями
элементов основных матриц, но и, что совсем необычно и непонятно, числом переменных состояния.
Из теории
следует, что число переменных состояния должно быть равно порядку выбранной системы дифференциальных
уравнений. Поэтому в системе (10), имеющей четвертый порядок, должно быть четыре переменных состояния.
В последнем случае, как видим, число переменных состояния возросло до десяти.
Резюмируя, отметим, что к процедурам создания lti -моделей относятся :
- ss -
создает модель пространства состояния по заданным матрицам A, B, C, D уравнений
состояния системы;
- dss - создает аналогичную модель по описанию пространства состояния более общего вида, когда
уравнения переменных состояния не разрешены относительно производных;
- tf - создает модель по заданным передаточным функциям системы;
- zpk -
создает модель по заданным нулям, полюсам и коэффициентам передачи системы;
- filt -
создает модель по дискретным передаточным функциям, записанным в форме полиномов
от
;
1
−
z
- set - присваивает значения некоторым другим полям LTI-объекта (таким, как названия входов и
выходов, название системы и т.п.).
Указанные процедуры позволяют создавать как непрерывные модели , так и дискретные. В последнем случае к
числу входных параметров процедуры следует добавить в конце значение параметра Ts - шага дискретизации, а
вводимые значения коэффициентов уже должны задавать
параметры дискретных передаточных функций (для
функций
tf и zpk), либо матрицы конечно-разностных уравнений пространства состояния - при
использовании процедур
ss и dss. При использовании процедуры filt должны задаваться векторы
227
коэффициентов числителя и знаменателя дискретной передаточной функции, представленной в виде отношения
полиномов от
. Приведем несколько примеров:
1
−
z
kzv1 = tf([1 4], [1 2 100])
Transfer function:
s + 4
---------------
s^2 + 2 s + 100
kzv2 = tf([1 4], [1 2 100],0.01)
Transfer function:
z + 4
---------------
z^2 + 2 z + 100
Sampling time: 0.01
kzv3 = tf([1 4], [1 2 100],'Variable','z^-1')
Transfer function:
1 + 4 z^-1
---------------------
1 + 2 z^-1 + 100 z^-2
Sampling time: unspecified
kzv4 =filt([1 4], [1 2 100])
Transfer function:
1 + 4 z^-1
---------------------
1 + 2 z^-1 + 100 z^-2
Sampling time: unspecified
Как следует из примеров, процедура filt полностью аналогична процедуре tf с добавлением в конец
списка входных параметров записи
'Variable'
,
'z^-1'
.
Те же процедуры
ss, dss, tf и zpk применяются также для преобразования моделей из одной из
указанных выше форм в другую. Первая и вторая применяются для преобразования модели в пространство
состояния, третья - в передаточную функцию, а четвертая - в нули-полюсы-коэффициент передачи.
Модель, заданную как
непрерывная система, можно перевести в дискретную форму, воспользовавшись
процедурой
c2d
в соответствии со схемой:
sysd = c2d(sys, Ts, method).
Здесь sys - исходная непрерывная заданная модель, sysd - получаемый в результате работы процедуры
дискретный аналог исходной системы, Ts - задаваемое значение шага дискретизации, method - параметр,
определяющий метод дискретизации. Последний параметр может принимать одно из таких значений:
- 'zoh' -
соответствует применению экстраполятора нулевого порядка: внутри интервала
дискретизации сигналы аппроксимируются постоянной величиной, равной значению сигнала в начале
интервала дискретизации;
- 'foh' -
соответствует применению экстраполятора первого порядка: внутри интервала
дискретизации сигналы аппроксимируются отрезками прямых, проходящих через концы кривой
сигнала в интервале дискретизации;
- 'tustin' - билинейная аппроксимация Тастина внутри интервала дискретизации;
- 'prevarp' - та же аппроксимация Тастина с заданной частотой предыскривления;
- 'matched' - метод согласования нуля и полюса.
Ниже приведены примеры перевода введенного ранее непрерывного колебательного звена kzv1 в дискретные
звенья по разным методам:
KZVd1= c2d(kzv1,0.01)
Transfer function:
228
0.01008 z - 0.009687
---------------------
z^2 - 1.97 z + 0.9802
Sampling time: 0.01
KZVd2= c2d(kzv1,0.01,'zoh')
Transfer function:
0.01008 z - 0.009687
---------------------
z^2 - 1.97 z + 0.9802
Sampling time: 0.01
KZVd3= c2d(kzv1,0.01,'foh')
Transfer function:
0.005029 z^2 + 0.0002308 z - 0.004864
-------------------------------------
z^2 - 1.97 z + 0.9802
Sampling time: 0.01
KZVd4= c2d(kzv1,0.01,'tustin')
Transfer function:
0.005037 z^2 + 0.0001975 z - 0.00484
------------------------------------
z^2 - 1.97 z + 0.9802
Sampling time: 0.01
KZVd5= c2d(kzv1,0.01,'prewarp',50)
Transfer function:
0.005145 z^2 + 0.000206 z - 0.004939
------------------------------------
z^2 - 1.97 z + 0.9798
Sampling time: 0.01
KZVd6= c2d(kzv1,0.01,'matched')
Transfer function:
0.01009 z - 0.009696
---------------------
z^2 - 1.97 z + 0.9802
Sampling time: 0.01
Процедура
d2c осуществляет обратную операцию - переводит дискретную систему в непрерывную, например:
k1 = d2c(KZVd1)
Transfer function:
s + 4
---------------
s^2 + 2 s + 100
k2 = d2c(KZVd4,'tustin')
Transfer function:
229
s + 4
---------------
s^2 + 2 s + 100
Как можно убедиться, указанные операции являются взаимно-обратными.
Процедура
d2d
позволяет переопределить дискретную систему, либо меняя шаг дискретизации
sys1 = d2d(sys,Ts),
либо вводя групповые задержки Nd (целое, в количестве шагов дискретизации)
sys1 = d2d(sys,[],Nd).
Приведем примеры. Вначале изменим шаг дискретизации на Ts=0.1 для системы KZVd1:
kd1=d2d(KZVd1,0.1)
Transfer function:
0.09352 z - 0.06018
-----------------------
z^2 - 0.9854 z + 0.8187
Sampling time: 0.1 ,
а затем введем задержку по входу, равную 3 Ts. Получим:
kd2=d2d(kd1,[],3)
Transfer function:
0.09352 z - 0.06018
-----------------------------
z^5 - 0.9854 z^4 + 0.8187 z^3
Sampling time: 0.1
Для создания модели нужно предварительно либо привести уравнения всей системы к форме уравнений
пространства состояний, либо найти передаточные функции системы. В общем случае это довольно сложная и
громоздкая задача. В то же время реальные системы автоматического управления (САУ) состоят из
соединенных между собой отдельных блоков (динамических звеньев), уравнения поведения которых обычно
достаточно просты. Поэтому в практике проектирования САУ принято использовать структурные методы,
когда САУ задается как определенная схема соединения отдельных элементарных динамических звеньев, и
фактически проектируется одно или несколько из этих звеньев таким образом, чтобы обеспечить заданное
качество всей системы. В соответствии с этим в MatLAB предусмотрена возможность «набирать» программно
«схему» САУ
путем предварительного ввода моделей звеньев, составляющих САУ, и последующего
«соединения» этих звеньев в единую структуру. К процедурам, осуществляющих расчет характеристик
соединений отдельных звеньев, относятся:
- plus (minus) - осуществляет «параллельное соединение» указанных в обращении звеньев, т. е.
определяет характеристики модели системы из параллельно соединенных звеньев; особенностью
является то, что вызов этих процедур может быть осуществлен не только обычным путем - указания
имени процедуры и перечисления (в скобках после имени) идентификаторов соединяемых звеньев, -
но и простым указанием идентификаторов звеньев, которые
должны быть объединены, с
простановкой между ними знаков « + » (при суммировании выходных сигналов звеньев) или « - »
(при вычитании выходных сигналов);
- parallel - осуществляет ту же процедуру параллельного соединения звеньев; в отличие от
предыдущей процедуры может использоваться для многомерных систем и осуществления
параллельного соединения лишь по некоторым входам и выходам;
- mtimes -
(или знак « * » между именами звеньев) - осуществляет последовательное соединение
звеньев, имена которых указаны; применяется для одномерных систем;
- series - последовательное частичное соединение многомерных систем;
- feedback - такое соединение двух звеньев, когда второе указанное звено составляет цепь
отрицательной обратной связи для первого звена;
- append - формальное объединение не связанных между собой систем (добавление выходов и
входов второй системы к выходам и входам первой);
230
- connect -
установление соединений выходов и входов многомерной системы, созданной
предварительно формальным объединением процедурой
append
; схема соединений задается
матрицей Q соединений, указываемой как один из входных параметров процедуры
- inv - рассчитывает САУ, обратную указанной, т. е. такую, у которой выходы и входы поменены
местами;
- vertcat - производит так называемую вертикальную конкатенацию (сцепление) систем (звеньев),
т. е. такое их объединение, когда входы их становятся общими, а выходы остаются независимыми;
для такого объединения необходимо, чтобы число входов объединяемых систем было одинаковым;
тогда число входов в результирующей системе останется таким же, как и каждой из объединяемых
систем, а число
выходов равно сумме выходов объединяемых систем;
- horzcat -
осуществляет «горизонтальное сцепление» указанных систем, при котором выходы
становятся общими, а входы добавляются.
Проиллюстрируем применение некоторых из этих процедур. Создадим модель углового движения торпеды
вокруг вертикали в виде последовательно соединенных двух звеньев: апериодического звена,
характеризующего влияние момента внешних сил относительно вертикали на угловую скорость торпеды
Torsk = tf(25,[100 50])
Transfer function:
25
----------
100 s + 50
и интегрирующего звена, описывающего переход от угловой скорости к углу поворота торпеды вокруг
вертикали
SkUg = tf(1,[1 0])
Transfer function:
1
-
s
Последовательное соединение этих звеньев можно осуществить двумя способами - применением процедуры
series
Tor1= series(Torsk,SkUg)
Transfer function:
25
--------------
100 s^2 + 50 s
либо просто операцией «перемножения» моделей
Tor = Torsk*SkUg
Transfer function:
25
--------------
100 s^2 + 50 s
Теперь сформируем цепь управления, входом которой является угол рыскания торпеды, а выходом - момент,
накладываемый на торпеду со стороны ее рулей направления. Ее будем предполагать состоящей из двух
параллельно соединенных частей - части, управляемой гироскопом направления и представляющей собой
обычное усилительное (статическое) звено
GN = tf(2,1)
Transfer function:
2
и части, управляемой гиротахометром, которую можно представить как дифференциально-колебательное звено
GT = tf([100 0],[1 10 100] )
Transfer function:
100 s
----------------