Волков В.Л. Моделирование процессов и систем в приборостроении
Подождите немного. Документ загружается.
61
Пример. Расчет дисперсии стохастического процесса на выходе динами-
ческой системы. Исходные данные: Sg(w) - спектральная плотность входного
процесса динамической системы, W(jw) - передаточная функция динамической
системы:
0201,198,1
25,0
)(
24
2
++
+
=
ww
w
wSg
, .
21243
122
)(
)(
2
++
+
==
ss
s
sg
sx
W(s)
1. Дисперсия стохастического процесса на выходе динамической системы
находится по спектральной плотности выходного процесса Sx(w). Определим
спектральную плотность выходного процесса по формуле из [1], представлен-
ной на стр.18: Sx(w)=Sg(w)|W(jw)|
2
. При этом рассчитаем сначала |W(jw)|
2
:
42
2
2
9450441
4144
)(
ww
w
jwW
++
+
=
.
Учитывая |W(jw)|
2
, Sg(w), найдем A
n
(s), B
n
(s) где s=jw:
.6132)122)(5,0()(
;21,2124,6603,72303)21243)(201,1()(
2
23422
++=++=
++++=++++=
sssssB
sssssssssA
n
n
Реализация алгоритма острема в Excel для порядка выше 2, также не вы-
зывает трудностей. Рассмотрим возмость построения алгоритма для n>5:
(
)
=
+
+
(
+
+
)
=
+
+
+
+
+
+
+
+
(
)
=
+
(
+
)
=
=
+
+
+
.
В этом случае надо заполнить формулами рекуррентного вычисления ко-
эффициентов a
k
, b
k
на значительно больший диапазон (например, для n=6). Это
представлено в следующей скрин-копии.
62
Рисунок 4.4. Дисперсия выхода в алгоритме Острема (Excel)
При использовании прикладной программы Astrom.exe исходные данные
в файле S(w).dat имеют следующий вид:
4 {порядок системы}
0 2 13 6 {полином B
n
(s)}
3 30 72.03 66.24 21.21 {полином A
n
(s)}
Результаты вычислений дисперсии получены в файле Astr.dat
k=4
3.00000 30.00000 72.03000 66.24000 21.21000
0.00000 2.00000 13.00000 6.00000
a[4]=0.10000 b[4]=0.00000
30.00000 65.40600 66.24000 21.21000
2.00000 13.00000 6.00000
D= 0.0000E+0000
a[3]=0.45867 b[3]=0.03058
65.40600 56.51153 21.21000
13.00000 5.35144
D= 1.0193E-0003
a[2]=1.15739 b[2]=0.23004
56.51153 21.21000
5.35144
D= 2.3881E-0002
a[1]=2.66438 b[1]=0.25231
21.21000
D= 3.5827E-0002
63
Результат вычислений дисперсии стохастического процесса на выходе
динамической системы получен, как и при расчетах в Excel (Dx=0,036).
4.3. Расчет частоты спектра и периода дискретности процесса
Оценивание максимальной частоты спектра стационарного стохастиче-
ского процесса и допустимого периода дискретности из условия заданной точ-
ности моделирования [20].
Максимальная частота спектра процесса (частота дискретизации) w
max
при динамической работе системы (в данном случае формирующего фильтра)
определяется при условии ослабления процесса с максимальной частотой по
отношению к процессам проходящим без ослабления до уровня погрешностей
при численном решении задачи моделирования процессов:
|W(jw
max
)|
2
=δ|W(jw
0
)|
2
, (4.8)
где w
0
- частота, для которой АФЧХ аппроксимируется горизонтальной
асимптотой, δ – относительный уровень ослабления сигнала на частоте дис-
кретизации.
Пример. Пусть процесс задан спектральной плотностью
S
g
(w)=(w
2
+0,25)/(w
4
+1,98w
2
+1,0201). Передаточная функция формирующего
фильтра для такого процесса имеет вид W
g
(s)=(s+0,5)/(s
2
+2s+1,01).
Зададим δ на уровне 10
-3
. Определим модуль амплитудо-фазочастотной
характеристики (АФЧХ) как отношение модуля числителя к модулю знамена-
теля. Учитывая
S
g
(w)=W
g
(jw)W
g
(-jw)=|W
g
(jw)|
2
,
получим для квадрата АФЧХ:
194,198,0
14
245,0
0201,198,1
25,0
)(
24
2
24
2
2
++
+
=
++
+
=
ww
w
ww
w
|w|Wg
.
При w
0
=0 АФЧХ не зависит от w, поэтому |W(jw
0
)|
2
=0,245. Для расчета
w
max
из (4.1) получим соотношение:
3
2
max
4
max
2
max
10
194,198,0
14
−
=
++
+
ww
w
,
64
откуда получим уравнение относительно w
max
0999,0998,300098,0
2
max
4
max
=−− ww ,
или 039,101959,4079
2
max
4
max
=−− ww . Из этого уравнения
2
max
w =4079,84 или w
max
=63,87 р/сек. В соответствии с импульсной теоремой
([1] стр.12) период дискретности процесса определяется в виде
dt=1/(2F
m
)=π/w
max
=0,049 c.
В соответствии с рассмотренной методикой расчета составлена програм-
ма Frequenc (табл.1.4). Входные данные (коэффициенты передаточной функ-
ции) задаются в файле W
g
(s).dat. В данном случае:
1 {число вариантов}
1 2 0.001 {порядок числителя, порядок знаменателя, погрешность}
0.5 1 {коэффициенты числителя}
1.01 2 1 {коэффициенты знаменателя}
Результаты расчета с помощью программы Frequenc.exe формируются в
файле Wc-out.dat и имеют следующий вид:
Computing Max Frequence dynamic system (Wc) and time of diskretization (dt)
Variant 1. Computing Wc,dt
---------------------------------------------------
W(s)-Function of dynamic system
5.00E-0001 1.00E+0000
1.01E+0000 2.00E+0000 1.00E+0000
K^2=0.2450740
|Wа(jw)|^2-Qwadrat of module AFFX Wp(jw)
(Co+,...CmW^(2m))/(Do+,...DnW(2n))
1.00E+0000 4.00E+0000
1.00E+0000 1.94E+0000 9.80E-0001
Biqwadrat equation of Wc
Wc^4-4078.420*Wc^2-1019.080=0
Qwadrat of poluse equation
-2.50E-0001 0.00E+0000 j
4.08E+0003 0.00E+0000 j
Max Frequence dynamic system Wc= 6.386E+0001
Time of diskretization dt= 4.919E-0002
65
4.4. Расчет непрерывной и дискретной матричной модели
Пример. Задана математическая модель объекта наблюдения в виде пере-
даточной функции следующего вида
2,116,168,1185,4
165,0
)()(
)(
)(
234
++++
+
=
+
=
ssss
s
svsu
sh
sW
.
Детерминированное входное воздействие u(t) задано суммой ступенчатой
и гармонической функций. Стохастическая составляющая g(t) определяется
формирующим фильтром D с белый шумом возмущения v(t) приведенном ко
входу (интенсивность шума возмущения Q=0,7).
Модель ОН имеет: X - вектор состояния системы; h(t) - сигнал выхода.
DvBuAXX ++=
&
,
,
CX
h
=
где A, B, D – переходные матрицы состояния (матрицы коэффициентов);
С - матрица линейного преобразования размера mxn.
=
0010
0001
C
.
Измеряются два параметра объекта наблюдения: выходной сигнал и его
производная. Уравнение измерений удовлетворяет линейному соотношению:
)()()( tNtCXtY
+
=
,
где N(t) - вектор шумов,
)()]()([М tRtNtN
Т
=
. Измерения производятся с шу-
мами (интенсивности шумов r
11
=13, r
22
=0,3).
Дискретная модель
)()()()1( kkkk
TvGuFXX
+
+
=
+
.
Матрицы F, G, T зависят от A, B, D Период дискретности dt для вычис-
ления этих матриц выбирается исходя из полосы пропускания динамической
системы в соответствии с импульсной теоремой [21].
66
Преобразование передаточных функций в матричные модели производят
стандартными методами, например, методом нормальной формы Коши.
Определение матричной модели ОН методом Коши производится с по-
мощью электронной таблицы Excel.
Результаты расчета непрерывной матричной модели представлены на
скрин-копии электронной таблицы (рисунок 4.5).
Рисунок 4.5. Таблица для расчета матриц непрерывной модели
Для расчета дискретной модели требуется определить период дискретно-
сти dt. На основании импульсной теоремы dt=π/ω
m
, где ω
m
– максимальная
частота пропускания системы. Вычисление периода дискретности произво-
дится в электронной таблице Excel (см. рисунок 4.6).
Рисунок 4.6. Реализация инструмента “Подбор параметра” для расчета dt
Метод расчета основан на определении ω
m
, когда затухание амплитудно-
частотной характеристики ограничено на уровне δ от уровня при ω = ω
0
= 0.
|
(
)
|
=
|
(
)
|
.
67
При δ = 10
-3
это уравнение приобретает вид
ω
+0,0775ω
+0,62ω
−10
ω
−10
=0.
Метод расчета основан на использовании инструмента Excel “Подбор па-
раметра”. В ячейке G8 набирается формула биквадратного уравнения, где в
качестве неизвестного параметра x = (ω
m
)
2
используется изменяемая ячейка
$F$8.
Рисунок 4.7. − Инструмент “Подбор параметра”
В данном случае формула биквадратного уравнения четвертого порядка в
ячейке G8 имеет вид
=F5+G5*F8+H5*F8^2+I5+F8^3+J5*F8^4.
Предварительно в искомой ячейке задается достаточно большое число,
заведомо большее, чем (ω
m
)
2
. После проведения подбора получены результаты,
представленные на рисунке 4.8.
Рисунок 4.8. Результаты подбора параметра x = (ω
m
)
2
68
В выделенной ячейке получена остаточная ошибка расчета, а в ячейках
H8, I8 найденные значения ω
m
и dt.
Дискретная модель содержит матрицы F, G, T (
dt
A
I
F
⋅
Ψ
+
=
,
dt
D
G
⋅
Ψ
=
).
Алгоритм расчета основан на представлении матричной экспоненты в виде
ряда Тейлора. При этом используется вспомогательная матрица
)!1(
,...,
2 +
+++=Ψ
k
dt
A
dt
AI
k
k
.
Итерационный расчет ψ производится до достижения условия точности
δ<
Ψ
Ψ−Ψ
+
+
)1(
)()1(
abs
k
kk
,
где δ - заданная относительная погрешность (δ=0,001);
||ψ
(k)
|| - норма матрицы ψ
(k)
, вычисляемая как сумма элементов матрицы.
Матрицы F, G, T, найденные с помощью электронной таблицы Excel,
представлены на копии экрана электронной таблицы (рисунок 4.9).
Рисунок 4.9. Результаты расчета матриц дискретной модели
69
5. Статистическое моделирование систем
Метод статистического моделирования, или метод Монте-Карло (МК), -
это численный метод решения математических задач при помощи моделиро-
вания случайных величин [22, 23]. Этот метод позволяет исследовать любой
процесс, на протекание которого влияют случайные факторы. Для задач, не
связанных со случайностями, можно искусственно придумать вероятностные
модели.
Общая схема применения метода состоит в том, что функцию, которую
нужно вычислить, представляют в виде математического ожидания f(x
1
,...,x
n
)
от n независимых случайных величин. Используя n случайных чисел, модели-
руемых датчиками случайных чисел, получают реализацию функции. Повто-
ряя, получают n случайных реализаций. В качестве приближенного значения
оцениваемой функции принимают среднее арифметическое по реализациям.
Среднеквадратическая погрешность полученной статистической оценки с рос-
том числа реализаций N убывает пропорционально 1/N
2
.
При статистическом моделировании систем производится воспроизведе-
ние процессов, происходящих в системах, с искусственной имитацией случай-
ных величин, от которых зависят эти процессы, с помощью датчиков случай-
ных чисел. Комбинируя детерминированные и стохастические зависимости,
составляют алгоритм моделирования системы. Применяя его, получают неза-
висимые реализации процесса в заданных условиях использования системы.
Характеристики, которые нужно определить, оцениваются методом МК. Алго-
ритм моделирования используется как для анализа, так и для оптимизации и
синтеза систем.
5.1. Оптимизация методом статистического моделирования
В качестве метода оптимизации функционала Q(F), зависящего от пара-
метров системы, используется метод статистического моделирования, легко
реализуемый на компьютере. При использовании этого метода, вектор пара-
метров F варьируют в допустимых границах случайным образом и вычисляют
при каждом статистическом испытании функционал Q(F). Затем при достиже-
нии достаточного статистического материала определяют экстремум функ-
ционала и выявляют соответствующие ему оптимальные значения вектора F.
Функционал часто определяется неявно при моделировании работы рас-
сматриваемой системы. Задача оптимизации решается при наличии ряда огра-
ничений на другие характеристики (функционалы Q(F)).
70
На рисунке 5.1 представлена структурная схема алгоритма оптимизации.
Рисунок 5.1. Структурная схема алгоритма оптимизации функционала
На схеме обозначено:
ДСЧ – датчик случайных чисел;
N – количество статистических испытаний;
D
P
– область допустимых второстепенных функционалов;
Q
min
(F) – минимальное значение ортимизируемого функционала;
F
опт
- вектор оптимальных параметров.
В качестве примера функционала можно, например, привести передаточ-
ную функцию ДПИ, зависящую от множества параметров, определяющих кон-
струкцию датчика и его динамические свойства. Большинство параметров
функционала являются неизменяемыми, но некоторые можно выбрать опти-
мальным образом для обеспечения экстремальных характеристик датчика, на-
-
-
k
=1
Вариация параметров
f
i
Расчет функционала
Расчет ограничений
P
(
F
)
P(F)∈D
Q
k
(
F
)<
Q
min
(
F
)
Q
min
(
F
)=
Q
k
(
F
)
k=k+1
-
k>N
ДСЧ
Q
min
(
F
) →
71
пример, по точности или по быстродействию. Изменяемыми параметрами век-
торного аргумента F в этом случае являются параметры корректирующего
устройства, выполняемого, например, на малогабаритном операционном уси-
лителе с навесными элементами.
Рассмотрим в качестве основного показателя качества динамики датчика
время переходного процесса t
ПП
, определяемое при воздействии скачкообраз-
ного сигнала на входе. Минимум времени переходного процесса, в данном
случае, обеспечивает минимум потерь полезной информации, т.е. минимум
интегрального показателя качества – среднеквадратической ошибки переход-
ного процесса. Другие важные показатели динамики (перерегулирование, сте-
пень затухания, запас устойчивости по фазе и амплитуде, полоса пропускания)
должны находиться в допустимых пределах, т.е. они являются функциями ог-
раничения. Таким образом, задача оптимизации (минимизации) времени пере-
ходного процесса поставлена как задача с ограничениями.
Показатель качества время переходного процесса определяется неявно,
т.е. при моделировании стохастической системы. В качестве системы рассмат-
рим ДПИ. Для моделирования ДПИ используется стандартная матричная мо-
дель, характеризующая динамическое движение. Модель определяется урав-
нением состояния с вектором состояния X=(x
1
,x
2
,…,x
n
) при входных воздейст-
виях U(t) и возмущениях V(t):
).()(
),()()()(
tCXty
tDVtButAXtX
=
++=
•
(5.1)
где A, B, C – матрицы коэффициентов соответствующей размерности.
Начальное состояние процессов в датчике линейных ускорений заданно
вектором X(0)=(x
10
,x
20
,…,x
n0
).
При поиске оптимальных параметров имеются ограничения, наклады-
ваемые на пределы изменения исходных независимых параметров и оптими-
зируемых выходных характеристик. Параметры F системы и второстепенные
показатели качества P
k
имеют ограничения:
∈
,
(
)
∈
, (5.2)
Моделирование динамики движения системы выполняется численным
решением матричных уравнений динамики (формула 5.1). Для этого состав-
ляются соответствующие дискретные матричные уравнения:
).()(
),()()1(
kCXky
kGukHXkX
=
+
=
+
(5.3)
72
где H, G – матрицы коэффициентов соответствующей размерности.
Большинство параметров системы, определяют ее конструкцию и стати-
ческую характеристику и являются неизменяемыми. Часть параметров, отно-
сящихся к изменяемой части (корректирующему устройству) может быть из-
менена в ограниченных пределах (определяемых в основном массо-
габаритными и конструктивными характеристиками).
5.2. Объектно-ориентированное статистическое моделирование
Решение примера оптимизации ДПИ проводится с использованием объ-
ектно-ориентированного программирования в среде Delphi [24]. С этой целью
разработана операционная оболочка программы, которая предназначена для
организации диалоговой настройки вычислительного эксперимента и его мно-
гократного проведения [25]. Основная панель (Desck Top) операционной обо-
лочки с кнопками управления представлена на рисунке 5.2.
Рисунок 5.2. Управляющая панель операционной оболочки
73
На основной панели размещены элементы управления:
1) Кнопка управления для запуска подпрограммы статистического моде-
лирования и построения графика процесса движения динамической системы.
2) Диалоговое окно для ввода и модификации исходных данных (незави-
симых переменных) исследуемой динамической системы. В этом окне преду-
смотрена возможность установки/сброса флажков разрешения случайной ва-
риации независимых переменных системы.
3) Движки, задающие: 1) число статистических испытаний; 2) макси-
мальное время наблюдения динамического движения системы.
4) Диалоговое окно для ввода и модификации векторной функции моде-
ли (зависимых переменных) от векторного аргумента (независимых перемен-
ных) исследуемой динамической системы.
5) Редактируемые окна для модификации порядка дифференциального
уравнения и размерности вектора независимых переменных исследуемой ди-
намической системы.
6) Диалоговое окно для просмотра результатов проектирования и опти-
мизации исследуемой динамической системы.
Кнопки управления. Кнопка “Моделирование” предназначена для запус-
ка подпрограммы статистического моделирования. При нажатии кнопки про-
изводится подготовка, проведение и завершение вычислительного экспери-
мента.
На этапе подготовки читаются исходные данные из входных файлов: по-
рядок полинома передаточной функции для числителя и знаменателя, число
исходных параметров системы, номинальные значения исходных параметров,
индикаторы варьирования переменных.
На этапе проведения вычислительного эксперимента производится мно-
гократный вычислительный эксперимент с выполнением в каждом статисти-
ческом испытании:
1) случайного варьирования исходных параметров системы x[1], x[2],...,
x[Nq] относительно номинальных значений в заданных диапазонах;
2) расчет коэффициентов b
0
, b
1
, b
2
, …, b
m
; a
0
, a
1
, …, a
n
передаточной
функции системы при исходных параметрах x[1], x[2],..., x[Nq];
3) расчет непрерывной матричной модели системы (формула 5.1);
4) расчет дискретной матричной модели системы (формула 5.3);
5) численное решение дифференциального уравнения системы при воз-
действии скачкообразного сигнала в заданном диапазоне времени;
6) определение времени переходного процесса, перерегулирования и
статистических характеристик этих показателей качества.
74
На этапе завершения вычислительного эксперимента (после проведения
заданного числа статистических испытаний) определяется минимальное время
переходного процесса (максимальное количество статистических испытаний в
данном пакете программ ограничено числом 1000) и соответствующие пара-
метры системы (оптимальные параметры). Строится график оптимального
переходного процесса в графическом окне (красная линия) и график переход-
ного процесса при номинальных значениях исходных параметров. Полученные
оптимальные параметры системы записываются в результирующий файл.
Кнопка “Результаты” позволяет вызывать диалоговое окно для просмотра
результатов проектирования и оптимизации исследуемой динамической сис-
темы (рисунок 5.3).
Рисунок 5.3. Диалоговое окно результатов проектирования
Движки числа статистических испытаний и времени наблюдения.
В графическом окне расположен Движок, устанавливающий максималь-
ное время наблюдения динамического движения системы (от 0 до 1,5 с).
Движок, задающий количество статистических испытаний системы рас-
считан на диапазон от 10 до 1000. Наиболее часто используемые значения (от
30 до 800) указаны синей полосой.
Диалоговое окно для просмотра результатов. В данном окне представле-
ны рассчитанные коэффициенты b
0
, b
1,
a
0
, a
1
, a
2
, a
3
, a
4
модели системы “Пере-
даточная функция”, результаты расчета матричной модели: матрицы A, B, C,
результаты проведенной оптимизации: минимальное время переходного про-
цесса и соответствующее ему перерегулирование и оптимальные значения па-
75
раметров ДЛУ. Результаты этого окна связаны с данными, записанными в вы-
ходные файлы: Model-W-Nominal.dat, Model-W-Optimal.dat, MK-0.dat, Koshi-
Model.dat. Также в файлах представлены дискретные модели системы для но-
минальных (Digital-Model-Nominal.dat) и оптимальных параметров (Digital-
Model-Optimal.dat). Эти файлы обновляются при моделировании.
Блок статистичекого моделирования, структурная схема алгоритма кото-
рого представлена на рисунке 5.4, содержит управляющую процедуру и блоки.
Рисунок 5.4. Блок статистичекого моделирования
76
Датчик случайных чисел. Случайные числа в соответствии с исходной
информацией о распределении аргументов вырабатываются программно с по-
мощью датчика случайных чисел (ДСЧ). При статистическом моделировании
ДС используются случайные величины, имеющие равномерное распределение.
Блок решения дифференциального уравнения состояния ДС. Этот блок
представлен процедурой Eyler(n,h,u,A,B,X), выполняющей численное решение
с высокой точностью. Блок имеет входные параметры: n – порядок динамиче-
ской системы (ДС), h – период дискретности по времени, u – установившееся
значение входного сигнала, A, B – матрицы коэффициентов решаемой системы
уравнений. Выходными параметрами являются компоненты вектора состояния
ДС: X.
Блок оптимизации показателя качества. Определяется минимальное вре-
мя переходного процесса, соответствующие ему перерегулирование и пара-
метры ДС.
Блок обработки результатов. Осуществляется накопление массивов,
оценку статистических характеристик, вычисление показателей качества, ре-
шение оптимизационных задач, графического представления результатов.
Процедуры вывода результатов проектирования и оптимизации. Про-
грамма реализует вывод данных в файлы (запись результирующих данных о
моделях ДС, времени переходного процесса и основных параметрах системы).
Это файлы: Model-W-Nominal.dat, Model-W-Optimal.dat, MK-0.dat, Koshi-
Model.dat. Также в файлах представлены дискретные модели системы для но-
минальных (Digital-Model-Nominal.dat) и оптимальных параметров (Digital-
Model-Optimal.dat). Эти файлы обновляются каждый раз при выполнении мо-
делирования.
5.3. Расчет передаточной функции ДС по заданным номинальным
параметрам и функциональным зависимостям
Расчет передаточной функции ДС выполняется программой Dynamic Sys-
tem в диалоговом режиме с использованием режима ввода/вывода с файлами.
При этом задаются номинальные значения параметров (до 50 параметров) и
функциональные зависимости коэффициентов передаточной функции от этих
параметров.
Вычисления производятся с помощью оригинальной программы, состав-
ленной на объектном Pascal. Вычисления производятся с использованием ре-
курсий и специальных функций. Для записи формулы используются круглые
скобки, арифметические знаки и специальные математические функции. Неза-
висимые переменные (аргументы) обозначаются как x[1],…, x[50].
77
Диалоговое окно для ввода и модификации исходных данных. В этом ок-
не производится изменение номиналов исходных параметров (независимых
переменных) исследуемой динамической системы. Здесь же с помощью
флажков-индикаторов (нулей и единиц) выбираются варьируемые параметры
(1 – разрешено варьирование, 0 – запрещено варьирование).
Исходные значения параметров и их индикаторные функции записаны в
файл Input.dat, который используется при вводе/выводе данных в программу.
Диалоговое окно для ввода и модификации векторной функции. В ниж-
нем диалоговом окне индицируются функциональные зависимости коэффици-
ентов передаточной функции: b
0
, b
1
, …, b
m
; a
0
, a
1
, …, a
n
(зависимых перемен-
ных) от исходных параметров: x[1], x[2],..., x[Nq] (независимых переменных)
исследуемой динамической системы.
Для ввода в программу функциональных зависимостей выражения запи-
сываются с соблюдением определенных правил в файле Function.dat. Правила
записи функциональных зависимостей индицируются при вхождении в диало-
говое окно для модификации этих зависимостей.
Окна для модификации порядка дифференциального уравнения и
размерности вектора независимых переменных. Передаточная функция ис-
следуемой динамической системы записана в общем виде (см. рисунок
5.2), где индексы m и n задаются дополнительно в редактируемых окнах.
Таким же образом указывается число исходных параметров Nq.
5.4. Преобразование модели “Передаточная функция”
в непрерывную и дискретную матричные модели
Для преобразования модели “Передаточная функция” в непрерывную и
дискретную матричные модели необходимо записать на основе передаточной
функции дифференциальное уравнение и затем преобразовать это уравнение в
систему дифференциальных уравнений 1-го порядка.
78
Динамическое движение системы характеризуется вектором состояния
X=(x
1
,x
2
,…,x
n
) при входном сигнале u(t) и уравнением выхода y(t), определяю-
щим выходной сигнал (формула 5.1).
Моделирование динамики движения системы выполняется численным
решением матричных уравнений динамики (формулы 5.1).
5.5. Численное решение дифференциальных уравнений
Моделирование динамики движения ДС выполняется численным реше-
нием системы уравнений динамики (формулы 5.3) методом Эйлера. Решение
относительно компонент x
i
выполняется с заданной точностью путем выбора
достаточно малого периода дискретности h:
'
(1)()()
iii
xkxkhxk
+=+
, (5.4)
где x’
i
- первая производная компоненты x
i
.
Для решения системы дискретных уравнений можно использовать элек-
тронную таблицу Excel, но для получения достаточно точного решения требу-
ется выбрать период дискретности на уровне 10
-6
секунды. С учетом того, что
время наблюдения движения ДС составляет от 0,001 с до 1 с, для получения
численных значений траектории движения ДС придется сформировать табли-
цу на уровне 1000000 значений. Решение такой задачи в Excel нерационально,
поэтому она решается при использовании Object Pascal и визуального про-
граммирования Delphi.
Для решения дифференциального уравнения, заданного моделью в виде
передаточной функции необходимо:
1) ввести коэффициенты передаточной функции;
2) сформировать непрерывную матричную модель уравнения состояния
(вычислить матрицы коэффициентов A, B, C);
3) задать период дискретности;
4) сформировать дискретную матричную модель уравнения состояния
(вычислить матрицы коэффициентов F, G);
5) используя процедуру пошагового решения системы разностных дис-
кретных уравнений, вычислить траекторию движения ДС с шагом дискретно-
сти по времени h (сформировать таблицу);
6) построить график движения ДС.
79
Решение дифференциального уравнения состояния ДС производится в
матричной форме с помощью процедуры:
Eyler(n:Integer;h,u:Real;A:Ms;B:Ms1;Var X:Ms1).
Входные данные процедуры:
n – порядок системы;
h – период дискретности;
u – входная переменная;
A, B – матрицы коэффициентов уравнения состояния.
Выходные данные процедуры:
X – вектор выходных переменных (вектор состояния).
Для примера рассматривается дифференциальное уравнение 4-го поряд-
ка, что не нарушает общности решения для более высокого порядка.
Дифференциальное уравнение имеет вид
)1(
100
)1(
1
)2(
2
)3(
3
)4(
4
)()( ubtubtyayayayaya +=++++
. (5.5)
Программа Dynamic System для данного примера дает диалоговое окно с
параметрами ДС и траекторию движения ДС (аналогично рисунку 5.2).
Все параметры моделей просматриваются в диалоговом окне “Результа-
ты” (аналогично рисунку 5.3) и фиксируются в файлах.
5.6. Метод статистического моделирования при определении
допусков на конструктивные параметры системы
Алгоритм исследования динамических систем методом МК состоит в
следующем: задаются случайные вариации исходным параметрам системы, по
ним находится значение целевой функции, процесс циклически повторяется с
накоплением массива, который по достижении нужного объема подвергается
статистической обработке.
Задача формирования допусков на исходные параметры ДС. Необходимо
обеспечить определенный уровень некоторого выходного свойства J ДС. Это
свойство является функцией n независимых аргументов F=(f
1
,…,f
n
)
T
J = ( f
1
,…, f
n
). (5.6)
На стадии проектирования ДС для J устанавливают его номинальное зна-
чение J
ном
и предельные допустимые отклонения – нижнюю J
1
и верхнюю J
2
80
границы. Они определяют полуширину поля допуска dJ=(J
2
-J
1
)/2 (в этом слу-
чае требования к уровню J задаются в виде
djJ ±
íîì
). Необходимо, чтобы J
попадало в поле допуска, тогда ДС по данному свойству считается работоспо-
собным. На основе полученных требований к параметру J решается обратная
задача − выявление допустимых значений аргументов, влияющих на J. Это
особенно важно для исходных параметров, имеющих высокие коэффициенты
влияния.
Установление численных значений границ J
1
и J
2
проводится на основе
тактико-технических требований к ДС. При статистическом моделировании
ДС находят среднеквадратическое отклонение и из условия dJ=k σ
j
определя-
ется σ
j
. (часто k - коэффициент ширины поля допуска принимают равным 3 из
условия доверительной вероятности 0,997).
Методом статистического моделирования находят коэффициенты чувст-
вительности (коэффициенты влияния) параметров f
i
ДС на выходную характе-
ристику. Определяют среднеквадратические отклонения аргументов и их до-
пустимые границы. Кроме поставленной задачи, вероятностной по своей при-
роде, более важной представляется математическая задача оптимизации. Здесь
случайные вариации исходных параметров вводятся искусственно. Задача оп-
тимизации предшествует задаче формирования допусков на параметры ДС.
Обеспечение допустимых отклонений характеристик ДС. Задача опреде-
ления коэффициентов чувствительности dJ/df
i
параметров ДС решается сле-
дующим образом. Производится моделирование работы ДС при номинальных
значениях параметров. При этом исследуемому параметру задаются случайные
вариации с заданным законом распределения. Регистрируется выходная харак-
теристика J и определяется ее дисперсия. Коэффициент чувствительности i-го
параметра определяется в виде
ifi
J
i
f
J
df
dJ
σ
σ
=
. (5.7)
Математическое обеспечение. В состав математического обеспечения
входят: 1) программа – модель ДС; 2) датчики случайных чисел; 3) программы
обработки результатов.
Моделирование ДС выполняется численным решением уравнения дина-
мики (5.1) методом прогноза и коррекции при начальном старте "грубым" ме-
тодом Эйлера.
Алгоритм вычислений. Алгоритм легко программируется для решения в
численном виде на ПЭВМ. Прикладная программа составлена в объектно-