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

Размещено на http://www.allbest.ru/

Введение

Статистическая динамика в конечном итоге определяет статистические характеристики выходных сигналов системы, с помощью которых можно оценить ее точность, надежность, пропускную способность. Это важные, но частные характеристики системы управления, ее эффективности.

Статистическая динамика и теория эффективности систем управления основываются на теории вероятностей, математической статистике и динамике систем.

Данная курсовая работа посвящена изучение «Статистического анализа непрерывных динамических систем». Он включает в себя три основных метода:

- метод эквивалентных возмущений;

- метод интерполяционных полиномов;

- метод статистических испытаний и его модификация.

Подробно рассмотрен метод интерполяционных полиномов. Все расчеты были проведены на персональном компьютере с помощью Borland Delphi.



1. Теоретическая часть

1.1 Методы статистического анализа нелинейных динамических систем

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

1.1.1 Метод эквивалентных возмущений

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

(1)

(2)

где- - n - мерный вектор неслучайных начальных условий; - к - мерный вектор центрированных некоррелированных между собой случайных величин, для которых известны моментные характеристики до q - порядка включительно,

т.е. (3)

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

(4)



В основе метода лежит гипотеза об известном характере зависимости между входными случайными параметрами и выходными величинами. Для рассматриваемого случая - это возможность ее представления в виде ряда Тейлора по случайным параметрам V1,… Vm до q - порядка включительно.

Пусть решения системы уравнений (1) являются функциями времени t и случайных величин V1,…,Vm :

(5)

Тогда рассматриваемая степень решения (5) может быть представлена в виде:

(6)

где (7)

Разложим функцию в ряд Тейлора по случайным параметрам V1,…,Vm в окрестности их нулевых математических ожиданий, ограничиваясь членами q -й степени:

+

(8)

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

- частные производные от функции по случайным величинам, вычисленные при их нулевых значениях. Применив (8) операцию отыскания математического ожидания, получим искомую вероятностную характеристику .

(9)

Однако решение (5) обычно неизвестно, поэтому найти частные производные в (9) затруднительно. Для устранения этой проблемы выбирается N различных комбинаций значений для входных случайных величин, называемых в дальнейшем эквивалентными возмущениями:

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

(10)

Сравнивая выражения (9) и (10), можно сделать вывод, что сумма Q будет равна , если выполнены условия:

; 

(11)

где

В результате решения системы уравнений (11) находятся неизвестные .

Исходная система уравнений (1) интегрируется N раз для найденных комбинаций эквивалентных возмущений, а затем рассчитываются и находятся искомые вероятностные характеристики по соотношениям вида:

(12)

Найдем решение системы (11) для случая q=2 (соотношения метода эквивалентных возмущений для q = 3 и q = 5).

Пусть для системы случайных величин V1,…,Vm известны вероятностные характеристики:

(13)

Тогда система (11) примет вид:

(14)

Для выбранного числа комбинаций и принятой матрицы эквивалентных возмущений из решения системы следует:

(15)

Таким образом, главные вероятностные характеристики выхода динамической системы M[xi] и D[xi] могут быть найдены по соотношениям:

(16)

(17)

Аналогичные формулы можно получить и при числе членов разложения в ряд Тейлора q > 2.

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

1).Выбор значений эквивалентных возмущений.

2).Численное решение системы дифференциальных уравнений (1) (m+2) раза для комбинаций значений случайных параметров, равных соответствующим эквивалентным возмущениям.

3).Вычисление искомых вероятностных характеристик по формулам типа (16), (17).

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

1.1.2 Метод интерполяционных полиномов

В основе интерполяционного метода лежит идея аппроксимации интерполяционными полиномами существующей зависимости между входными случайными параметрами и выходными величинами. Метод также применяется при математическом описании динамической системы в виде (1). При этом считается, что случайные величины V1…Vm являются независимыми и для них известны законы распределения p(V1)… p(Vm).

Пусть решение системы (1) имеет вид:

(1)

Требуется найти оценку математического ожидания некоторой функции от фазовых координат системы в момент времени t (под эту формулировку подходят любые вероятностные характеристики):

= (2)

Так как явный вид выражения (1) неизвестен, то можно найти значения решений (1) путем численного интегрирования (1) на ЭВМ при определенных значениях случайных параметров V1…Vm, которые называются узлами интерполирования.

Пусть для первой случайной величины V1 выбраны определенным образом q1 вариантов узлов интерполирования ; для произвольной j -й случайной величины Vj выбраны qj вариантов узлов интерполирования ; для последней m -й случайной величины Vm выбраны qm вариантов узлов интерполирования .

Тогда общее число расчетных вариантов .

Представим функцию приближенно с помощью интерполяционного полинома Лагранжа:

(3)

где - многочлен степени qj относительно переменной Vj,

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

Тогда искомую характеристику (3) М[Ф(...)] можно представить с учетом (2), поменяв местами порядок интегрирования и суммирования, в виде:

(4)

Здесь:

(5)

Называются числами Кристоффеля.

Доказано, что если в качестве узлов интерполирования выбрать корни ортогональных полиномов по весу, равному плотности распределения случайной величины p(V), то при использовании n узлов интерполирования данный метод дает точные значения в классе многочленов всех степеней до степени q=2n-1 включительно. При этом для каждой плотности распределения вероятностей, являющейся весовой функцией, существует единственная система ортогональных многочленов.

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

Пусть случайная величина V распределена равномерно на интервале [a,b]. Тогда с помощью линейного преобразования вида:

(6)

случайную величину V можно получить из случайной величины , распределенной равномерно на интервале [-1, 1] с плотностью .

Ортогональной системой многочленов на интервале [-1, 1] с весом Ѕ является система многочленов Лежандра вида:

(7)

Числа Кристоффеля для узлов, выбранных в точках, соответствующих корням многочленов Лежандра степени (n+1), определяются по формуле:

(8)

Значения корней ортогональных многочленов Лежандра до 8 порядка включительно (q<=8) для случая закона равномерного распределения плотности вероятности на интервале [-1, 1]. Их следует пересчитывать для конкретных параметров равномерного закона распределения случайной величины по (6).

(9)

Числа Кристоффеля pk рассчитываются заранее и для q<=8, при этом они выбираются из таблицы без пересчета.

Пусть случайная величина V распределена по нормальному закону с параметрами MV, DV. Тогда с помощью линейного преобразования вида:

(10)

случайную величину V можно получить из нормально распределенной случайной величины с, .

Весовой функции видасоответствует система ортогональных многочленов Эрмита на интервале , которые можно вычислить по рекуррентным формулам:

(11)

Значения корней ортогональных многочленов Эрмита для веса pr() и пересчитываются для случайной величины V по (10):

(12)

Числа Кристоффеля pr не зависят от параметров нормального закона распределения. Данные могут быть рассчитаны и для других законов распределения случайных величин.

Опишем этапы статистического анализа динамической системы при использовании интерполяционного метода:

1.Для заданных порядков интерполяционных многочленов q1,…,qm по случайным величинам V1,…,Vm узлы интерполирования и числа Кристоффеля .

2. По формулам пересчета типа (6), (12) рассчитываются значения узлов интерполирования для каждой случайной величины - , а также числа Кристоффеля вида (5).

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

4. Вычисляется искомая вероятностная характеристика по соотношению (4).

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

Если начинает выполняться неравенство

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

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

1.1.3 Метод статистических испытаний и его модификация

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

Сущность метода статистических исследований (МСИ) заключается в следующем. На первом этапе создается моделирующий алгоритм, имитирующий воздействие случайных и детерминированных факторов на процесс функционирования и определяющий интересующие характеристики динамических систем:

(1)

(2)

Далее многократно повторяются реализации алгоритма (выполняется моделирование динамической системы, при которых случайные факторы V, о, x0 заменяются их реализациями, полученными случайным образом) и формируется статистическая выборка (второй этап):

(3)

(4)

где - реализации случайных факторов V, о , x0 соответственно, N - объем выборки реализации.

Затем обрабатывается полученная статистическая выборка реализаций методами математической статистики, и находятся искомые вероятностные характеристики выхода динамической системы (третий этап):

(5)

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

Так как в (5) оценка является суммой случайных величин, равнозначных в вероятностном отношении, то согласно центральной предельной теореме Ляпунова закон распределения случайной величины будет при N>? приближаться сколь угодно близко к нормальному закону распределения. Поэтому в качестве меры относительной погрешности оценки можно принять величину:

(6)

Тогда можно записать следующую вероятностную оценку погрешности:

(7)

где ? - заданная верхняя граница погрешности дN, - вероятность не превышения погрешностью дN фиксированной границы ?, Ф(…) - функция Лапласа.

Из (7) следует, что если гарантийную вероятность принять равной Рг=0,99, то формула для определения необходимого числа опытов N для обеспечения требуемой точности e с вероятностью Рг запишется в виде:

откуда: (8)

Таким образом, определение оценки вероятностных характеристик по методу статистических испытаний требует проведения большого объема моделирований (например, достижение 5% -й точности требует проведения не менее 400 испытаний), что является главным недостатком упомянутого метода.

Для снижения трудоемкости метода статистических испытаний разработаны различные его модификации. Одна из наиболее простых и в то же время эффективных модификаций основана на дополнительной обработке выборки входных случайных факторов и случайных выходных контролируемых характеристик.

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

Если математическая модель динамической системы преобразована к виду (1), то уточненная оценка вероятностной характеристики M*[yi] определяется по соотношению:

(9)

Здесь M[V] - точное значение вектора математических ожиданий входных случайных величин V;

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

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



2. Практическая часть

2.1 Постановка задачи

Пусть имеется некоторое линейное множество R действительных функций, определенных на отрезке [a, b], и некоторая конечная или счетная совокупность достаточно простых функций этого множества {ji(x)}. Предполагается, что любая конечная подсистема функций {ji(x)} линейно независима на [a, b]. Рассматривается задача приближения функции f(x) R функциями вида Qm(x) = c0j0(x) + c1j1(x) +...+ cmjm(x). На [a, b] выбирается некоторая фиксированная совокупность попарно различных точек x0, x0,..., xn и каждой конкретной функции f(x) R ставится в соответствие обобщенный многочлен Qm(x), значения которого в выбранных точках совпадают со значениями функции f(x). Точки x0, x0,…, xn называются узлами интерполирования, а обобщенный многочлен, обладающий указанным свойством, называется обобщенным интерполяционным многочленом для f(x) по заданной системе узлов.

2.2 Алгоритм решения задачи

Задаётся начальная функция и интервал, в пределах которого строится данная функция. Затем следует построить заданную функцию по выбранным произвольно точкам. Далее по формуле полинома Лагранжа в выбранных узлах интерполирования, рассчитывается и строится, приближенная к заданной, функция.

На выходе программы должны получить два графика:

- сама функция;

- интерполяционный многочлен.

Код программы представлен в приложении А, результат работы программы представлен в приложении Б.

3. Техническая часть

3.1 Обоснование выбора программного обеспечения

Выбор среды разработки обусловлен тем, что посредством Delphi имеется возможность написания наглядной программы с дружественным интерфейсом и большим количеством возможностей. Синтаксис Borland Pascal и Borland Delphi очень похож, и среда разработки Borland Delphi становится всё популярнее у разработчиков, так как охватывает обширную часть методов и технологий для взаимодействия с операционной системой и оборудованием компьютера.

Delphi - это среда быстрой разработки, в которой в качестве языка программирования используется язык Delphi. Язык Delphi - строго типизированный объектно-ориентированный язык.



Вывод

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

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



Список используемой литературы

1. Болнокин В.Е., Чинаев П.И. «Анализ и синтез САУ на ЭВМ» стр.45;

2. Семенов В.В., Пантелеев А.В., Бортаковский А.С. «Математическая теория управления в примерах и задачах»;

3. Жданюк Б.Ф. «Основы статистической обработки траекторных измерений»;

4. Астапов Ю.М., Медведев В.С. «Статистическая теория систем автоматического регулирования и управления»;

5. http://www.ssau.ru/resources/ump/belokonov_sads/6/.



Приложение А

статистический нелинейный интерполяционный полином

(Листинг программы)

unit MainUnit;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls, ExtCtrls, ComCtrls, TeEngine, Series, TeeProcs, Chart;

type

Tf1 = class(TForm)

bpaint: TButton;

ledita: TLabeledEdit;

leditb: TLabeledEdit;

ComboBoxn: TComboBox;

ComboBoxfunction: TComboBox;

butassessment: TButton;

panel: TPanel;

Label1: TLabel;

Label2: TLabel;

im2: TImage;

Chart1: TChart;

Series1: TFastLineSeries;

Series2: TFastLineSeries;

procedure bpaintClick(Sender: TObject);

procedure ComboBoxnChange(Sender: TObject);

procedure FormCreate(Sender: TObject);

procedure butassessmentClick(Sender: TObject);

procedure ComboBoxfunctionChange(Sender: TObject);

procedure FormResize(Sender: TObject);

procedure FormShow(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

Приложение А

(Обязательное)

end;

var

f1: Tf1;

implementation

uses

ReadFileUnit;

{$R *.dfm}

type

elem = extended;

const

maxn = 20;

var

x, y: array[0..maxn + 1] of elem;

rr: array[1..maxn + 1, 1..maxn + 2] of elem;

xBegin, xEnd : Real;

n: byte;

function f(const x: elem): elem;

begin

if f1.ComboBoxfunction.ItemIndex = 0 then f := sin(pi * x)

else if f1.ComboBoxfunction.ItemIndex = 1 then f := cos(pi * x)

else f := (x * x - 1) * sin(4 * x);

end;

function Lagrange(const x0: elem): elem;

var

s, p: elem;

i, j: byte;

begin

(Обязательное)

s := 0;

for i := 0 to n do

begin

p := 1;

j := 0;

while j < i do

begin

p := p * (x0 - x[j]) / (x[i] - x[j]);

inc(j);

end;

j := i + 1;

while j <= n do

begin

p := p * (x0 - x[j]) / (x[i] - x[j]);

inc(j);

end;

s := s + y[i] * p;

end;

result := s;

end;

function Polynom(const x: elem): elem;

begin

result := Lagrange(x)

end;

procedure Tf1.bpaintClick(Sender: TObject);

var

i, n : Integer;

d, mn, mx, y : Real;

begin

f1.ComboBoxnChange(f1.ComboBoxn);

if ledita.Text = '' then ledita.Text := '0';

if (leditb.text = '') or (strtofloat(ledita.Text) >= strtofloat(leditb.Text)) then

leditb.Text := floattostr(strtofloat(ledita.Text) + 1);

Приложение А

(Обязательное)

Chart1.LeftAxis.Automatic := True;

Chart1.BottomAxis.Automatic := False;

Series1.Clear;

Series2.Clear;

Series1.SeriesColor := clGreen;

Series2.SeriesColor := clRed;

xBegin := StrToFloat(ledita.Text);

xEnd := StrToFloat(leditb.Text);

Chart1.BottomAxis.Minimum := 0.99*xBegin;

Chart1.BottomAxis.Maximum := 1.01*xEnd;

n := Chart1.Width;

if n > 100 then n := 100;

d := (xEnd-xBegin) / n;

mn := f(xBegin);

mx := mn;

while xBegin < (xEnd+d) do

begin

y := f(xBegin);

Series1.AddXY(xBegin, y, '', clGreen );

if y < mn then mn := y;

if y > mx then mx := y;

y := Polynom(xBegin);

Series2.AddXY(xBegin, y, '', clRed );

if y < mn then mn := y;

if y > mx then mx := y;

xBegin := xBegin + d;

end;

Chart1.LeftAxis.Automatic := False;

Chart1.LeftAxis.Minimum := 0.99*mn;

Приложение А

(Обязательное)

Chart1.LeftAxis.Maximum := 1.01*mx;

end;

procedure Tf1.ComboBoxnChange(Sender: TObject);

var

h: elem;

i, j: byte;

begin

n := strtoint(f1.comboboxn.Items[f1.ComboBoxn.ItemIndex]);

n := strtoint(comboboxn.Text);

if ledita.Text = '' then ledita.Text := '0';

if (leditb.text = '') or (strtofloat(ledita.Text) >= strtofloat(leditb.Text)) then

leditb.Text := floattostr(strtofloat(ledita.Text) + 1);

h := (strtofloat(leditb.Text) - strtofloat(ledita.Text)) / n;

for i := 0 to n do

begin

x[i] := strtofloat(ledita.Text) + i * h;

y[i] := f(x[i]);

end;

for i := 1 to n + 1 do

begin

rr[i, 1] := strtofloat(ledita.Text) + (i - 1) * h;

rr[i, 2] := f(rr[i, 1]);

end;

for j := 3 to n + 2 do

for i := 1 to n + 3 - j do rr[i][j] := (rr[i][j - 1] - rr[i + 1][j - 1]) / (rr[i][1] - rr[i + j - 2, 1]);

end;

procedure Tf1.FormCreate(Sender: TObject);

var

i: byte;

begin

Series1.Clear;

Series2.Clear;

f1.Caption := application.Title;

comboboxn.ItemIndex := 0;

f1.ComboBoxnChange(comboboxn);

comboboxn.Items.Clear;

for i := 2 to maxn do comboboxn.Items.Add(inttostr(i));

comboboxn.ItemIndex := 0;

bpaint.Click;

im2.Canvas.Pen.Color := clBlack;

im2.Canvas.Brush.Color := clGreen;

im2.Canvas.FillRect(Rect(5, 5, 13, 13));

im2.Canvas.Brush.Color := clWhite;

im2.Canvas.TextOut(15, 1, ' --- Функция f(x)');

im2.Canvas.Brush.Color := clRed;

im2.Canvas.FillRect(Rect(5, 25, 13, 33));

im2.Canvas.Brush.Color := clWhite;

im2.Canvas.TextOut(15, 21, ' --- Интерполяционный многочлен');

end;

procedure Tf1.butassessmentClick(Sender: TObject);

const

h = 1000;

var

max, xm, x, a, b, r: elem;

i,

ind: integer;

g: textfile;

s: string;

begin

s := 'Assessment_Lagrange.txt';

assignfile(g, s);

rewrite(g);

if f1.ledita.Text = '' then f1.ledita.Text := '0';

if (f1.leditb.text = '') or (strtofloat(f1.ledita.Text) >= strtofloat(f1.leditb.Text)) then

f1.leditb.Text := floattostr(strtofloat(f1.ledita.Text) + 1);

a := strtofloat(f1.ledita.Text);

b := strtofloat(f1.leditb.Text);

xm := (b - a) / h;

writeln(g, 'f(x)=', f1.ComboBoxfunction.Text, ' x belong [', a:7, ',', b:7, ']');

writeln(g, 'n':3, ' ', 'Assessment');

writeln(g);

ind := f1.ComboBoxn.ItemIndex;

f1.ComboBoxn.ItemIndex := 0;

f1.ComboBoxnChange(f1.ComboBoxn);

for i := 2 to maxn do

begin

x := a;

max := -1;

while x <= b do

begin

r := abs(f(x) - Polynom(x));

if r > max then max := r;

x := x + xm;

end;

writeln(g, n:3, ' ', max:11);

f1.ComboBoxn.ItemIndex := f1.ComboBoxn.ItemIndex + 1;

f1.ComboBoxnChange(f1.ComboBoxn);

end;

f1.ComboBoxn.ItemIndex := ind;

f1.ComboBoxnChange(f1.ComboBoxn);

closefile(g);

f2.Memo1.Lines.LoadFromFile(s);

f2.Show;

f2.Caption := s;

end;

procedure Tf1.ComboBoxfunctionChange(Sender: TObject);

begin

f1.ComboBoxnChange(comboboxn);

end;

function podgon(s: string): string;

begin

while length(s) < 25 do s := s + ' ';

podgon := s;

end;

procedure Tf1.FormResize(Sender: TObject);

begin

bpaint.Click;

end;

procedure Tf1.FormShow(Sender: TObject);

begin

ComboBoxfunction.ItemIndex := 0;

bpaint.Click;

end;

end.

unit ReadFileUnit;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls;

type

Tf2 = class(TForm)

Memo1: TMemo;

procedure FormShow(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

f2: Tf2;

implementation

{$R *.dfm}

procedure Tf2.FormShow(Sender: TObject);

begin

Memo1.Align := alClient;

end;

end.



Приложение Б

Рисунок 1 - Результат работы программы

Рисунок 2 - Расчеты по формуле Лагранжа

Размещено на Allbest.ru

...