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

1. Аналитическое решение

2. Стандартная схема статистического моделирования

3. Рациональная схема статистического моделирования Заключение Список использованных источников Приложения

Требуется определить математическое ожидание выходного сигнала X неустойчивого апериодического звена в заданный момент времени Т. Модель звена:

где g = G (t), X (0) = A.

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

Допустимая абсолютная погрешность результата: е_доп. = 0,01.

Задачу решить тремя способами:

· Используя стандартную схему статического моделирования;

· Используя рациональную схему статистического моделирования с применением метода расслоенной выборки;

· Аналитически.

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

При использовании рациональной схемы статистического моделирования обеспечить снижение требуемого количества опытов по сравнению со стандартной схемой не менее чем в 10 раз.

Исходные данные (вариант 2−2):

G = 1? 1.4,

a = 0.6? 0.8,

T = 1.3,

A = 1,

k = 1.2.

1. Аналитическое решение

статистический моделирование математический апериодический Решим дифференциальное уравнение вида[1]:

(1)

где g = G (t),

X (0) = A.

Сначала найдем решение соответствующего однородного дифференциального уравнения:

Подставим полученное решение однородного дифференциального уравнения в (1):

Найдем С₁ из условия X (0) = A:

В результате имеем:

Решение исходного дифференциального уравнения (1) имеет вид:

(2)

где g — случайный параметр, распределенный по равномерному закону в интервале [1;1.4],

a — случайный параметр, распределенный по равномерному закону в интервале [0.6;0.8],

Для Т=1.3 с учетом статистической независимости k и g определим искомую характеристику:

где — искомое математическое ожидание.

С учетом (1) находим :

Таким образом,

Определим дисперсию :

(3)

где — дисперсия выходного сигнала.

Введем обозначение: и найдем :

(4)

Рассчитаем слагаемые, входящие в (4):

;

Таким образом, 21.77.

Подставив полученные значения в (3), определим дисперсию выходного сигнала:

С учетом известной дисперсии оценим необходимое количество опытов с погрешностью :

где — необходимое количество опытов.

Значение параметра зависит от доверительной вероятности. Примем P_д=0,997 и a_д=3. Подставив значения параметров в (5), получим:

опытов.

Все перечисленные расчеты производились в математическом пакете MathCAD [2], приводятся в Приложении А.

2. Стандартная схема статистического моделирования

Если трудоемкость эксперимента имеет существенное значение, применяются итерационные алгоритмы получения оценок. Идея итерационных алгоритмов состоит в том, что определение точности и требуемого количества опытов проводится в ходе эксперимента на основе получаемых оценок искомых параметров. Блок-схема типового итерационного алгоритма приведена на рисунке 1.

Рисунок 1 — Блок-схема итерационного алгоритма

Для задачи оценки математического ожидания случайной величины x предусматривается:

1. Проведение начальной серии опытов объемом n и накопление сумм

где — реализация случайной величины x в отдельных опытах.

2. Вычисление оценок математического ожидания и дисперсии :

(6)

. (7)

3. Получение оценки требуемого количества опытов:

. (8)

4. Проведение дополнительной серии опытов объемом и накопление сумм:

.

5. Уточнение оценок математического ожидания m^*_x и дисперсии D^*_x:

(9)

. (10)

Провели начальную серию опытов n = 200. Накопили суммы и: Вычислили оценки математического ожидания и дисперсии по (6) и (7): Получили оценку требуемого количества опытов по (8): Так как, то провели дополнительную серию опытов Для того, чтобы не проводилось лишнее число опытов искусственно уменьшили n' в 2 раза. Таким образом, опытов. Вновь накопили суммы, и уточнили оценки математического ожидания и дисперсии по (9) и (10): Тогда оценка требуемого количества опытов получилась: Значение n = 16 260+200=16 460 опытов.

После данной итерации 16 460<22 806, следовательно, продолжили выполнение итерационного алгоритма. Получили следующие результаты:

.

Проверили выполнение условия. Данное условие не выполнилось, так как 22 806>22 685, следовательно, алгоритм завершил работу.

Окончательные результаты :

Дифференциальное уравнение (1) решается численным интегрированием методом Эйлера первого порядка с шагом 0.001. Программа, реализующая итерационный алгоритм, написана в среде Borland Delphi 7. Текст программы представлен в Приложении Б.

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

· Выборочный закон распределения может существенно отличаться от нормального. Чаще всего оценки требуемого количества опытов оказываются завышенными.

· Разброс составляющих выборку реализаций случайной величины может оказаться существенно меньше истинного ее разброса.

· Оценки требуемого количества опытов оказываются резко заниженными, а результаты моделирования — неточными. Во избежание подобных ситуаций рекомендуется выбирать объем начальной серии опытов не менее 100−500.

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

3. Рациональная схема статистического моделирования

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

В соответствии с данным методом область G возможных значений случайного вектора разбивается на K=10 непересекающихся областей G_k:

Метод предполагает проведение статического моделирования для каждой из областей G_k с использованием для вектора случайных параметров плотностей распределения вероятностей где p_k — вероятность попадания случайного вектора V в область G_k

.

В нашем случае p_k = 0.1.

Блок-схема итерационного алгоритма метода расслоенной выборки приведена на рисунке 2.

Рисунок 2 — Блок-схема итерационного алгоритмаметода расслоенной выборки.

1.Если для области G_k выполним N_k опытов, получим оценку математического ожидания искомого показателя для данной области:

. (11)

Результирующая оценка должна рассматриваться как случайная дискретная величина, значения которой наблюдаются с вероятностями p_k. Тогда результирующая оценка определяется усреднением:

. (12)

2.Определим дисперсию оценки (9), имея в виду, что все N₁+ N₂+ N₃+…+ N₁₀ слагаемые — независимые случайные величины:

. (13)

Дисперсия случайной величины может быть оценена следующим образом:

(14)

3. Введя в рассмотрение доли от общего количества опытов, соответствующие областям G_k,

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

(15)

При удачном разбиении области G и удачном выборе соотношения количества опытов для отдельных областей G_k дисперсия оценки (13) может быть существенно снижена. Оптимальные значения должны быть пропорциональны произведениям

Провели начальную серию опытов N=200. После проведения данной серии опытов были получены следующие результаты:

· Оценка математического ожидания для каждой из 10 областей на основании (11):

· Результирующая оценка математического ожидания по (12):

· Дисперсия для каждой из 10 областей по (14):

· Дисперсия оценки математического ожидания по (13):

· Требуемое количество опытов, рассчитанное по (15):

опытов.

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

После второй итерации получили:

· N = 1982 опытов.

· Оценка математического ожидания для каждой из 10 областей:

· Результирующая оценка математического ожидания:

· Дисперсия для каждой из 10 областей:

· Дисперсия оценки математического ожидания:

(

· Требуемое количество опытов:

опытов.

После третьей итерации алгоритма:

· N = 2191 опытов.

· Оценка математического ожидания для каждой из 10 областей:

· Результирующая оценка математического ожидания:

· Дисперсия для каждой из 10 областей:

· Дисперсия оценки математического ожидания:

· Требуемое количество опытов:

опытов.

Дифференциальное уравнение (1) решается численным интегрированием методом Эйлера первого порядка с шагом 0.001. Программа, реализующая данный метод снижения трудоемкости, написана на языке Delphi 7.

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

Заключение

По заданию работы требовалось определить математическое ожидание выходного сигнала X апериодического звена в момент времени T тремя методами. В результате решения данной задачи тремя способами были получены следующие результаты:

§ Используя стандартную схему статистического моделирования

§ Используя метод расслоенной выборки

§ Аналитически

Использование метода расслоенной выборки обеспечило снижение требуемого количества опытов по сравнению со стандартной схемой в 10.85 раз.

Список использованных источников

1. Бертмант А. Ф. Краткий курс математического анализа. — М.: Наука, 1965.

2. Кирьянов Д. Самоучитель MathCAD 11. СПБ.: Бхв-Петербург, 2003.

3. Емельянов В. Ю. Методы моделирования стохастических систем управления. — СПб.: БГТУ, 2004.

4. Потапов М. К. Алгебра и анализ элементарных функций. — М.: Наука, 1980.

5. Бобровский С. DELPHI 7. Учебный курс.- СПБ.: Питер, 2003.

Приложения

Приложение А

Аналитические расчеты, произведенные в математическом пакете MathCAD.

Приложение Б

Программная реализация стандартной схемы статистического моделирования

unit Standart;

interface

uses

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

Dialogs, StdCtrls, Buttons;

type

TForm1 = class (TForm)

Memo1: TMemo;

Edit1: TEdit;

Memo2: TMemo;

Button1: TButton;

BitBtn1: TBitBtn;

procedure Button1Click (Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

const

alfa=9;

h=0.001;

eps=0.01;

kk=1300;

k=0.5;

var

Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1. Button1Click (Sender: TObject);

var

Sx, Sx2, mx, Dx, G, a, x, Sx2n, Sxn, mxn, Dxn: real;

i, j, ll:integer;

ntreb, nn, n, n1,n2:longint;

flag:boolean;

begin

Sx:=0;

Sx2:=0;

randomize;

nn:=200;

Memo1.Lines.Clear;

Memo2.Lines.Clear;

for i:=1 to nn do

begin

G:=random*0.4+1;

a:=random*0.2+0.6;

x:=1;

for j:=1 to kk do x:=x+a*x*h+g*1.2*a*h;

Sx:=Sx+x;

Sx2:=Sx2+sqr (x);

end;

mx:=Sx/nn;

Dx:=Sx2/(nn-1)-sqr (mx);

ntreb:=round (alfa*Dx/sqr (eps));

Memo1.Lines.Add ('Mx='+FloatToStr (mx));

Memo1.Lines.Add ('Dx='+FloatToStr (Dx));

Memo1.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo1.Lines.Add ('Sx='+FloatToStr (Sx));

Memo1.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo1.Lines.Add ('nn='+FloatToStr (nn));

n1:=0;

flag:=true;

while nn

begin

n1:=ntreb-nn;

edit1.text:=floattostr (n1);

if (n1>8500) then n1:=n1 div 2

else begin

if n1<7 then n1:=n1*2;end;

for i:=nn to n1+nn do

begin

G:=random*0.4+1;

a:=random*0.2+0.6;

x:=1;

for j:=1 to kk do x:=(a*x+g*1.2*a)*h+x;

Sx:=Sx+x;

Sx2:=Sx2+sqr (x);

end;

Memo2.Lines.Add ('n1='+FloatToStr (n1));

mx:=(Sx)/(n1+nn);

Dx:=(Sx2)/(n1+nn-1)-sqr (mx);

Memo2.Lines.Add ('Sxn='+FloatToStr (Sx));

Memo2.Lines.Add ('Sx2n='+FloatToStr (Sx2));

Memo2.Lines.Add ('mx='+FloatToStr (mx));

Memo2.Lines.Add ('Dx='+FloatToStr (Dx));

nn:=n1+nn;

ntreb:=round (alfa*Dx/sqr (eps));

Memo2.Lines.Add ('nn='+FloatToStr (nn));

Memo2.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo2.Lines.Add ('************************');

end.

Приложение В

Программная реализация метода расслоенной выборки

unit Viborka;

interface

uses

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

Dialogs, Buttons, StdCtrls;

type

TForm1 = class (TForm)

Memo1: TMemo;

Memo2: TMemo;

Memo3: TMemo;

Memo4: TMemo;

Memo5: TMemo;

Memo11: TMemo;

Memo12: TMemo;

Button1: TButton;

BitBtn1: TBitBtn;

Edit1: TEdit;

Edit2: TEdit;

Edit3: TEdit;

StaticText1: TStaticText;

StaticText2: TStaticText;

StaticText3: TStaticText;

StaticText4: TStaticText;

StaticText5: TStaticText;

StaticText6: TStaticText;

StaticText7: TStaticText;

StaticText8: TStaticText;

StaticText9: TStaticText;

StaticText10: TStaticText;

Memo6: TMemo;

Memo7: TMemo;

Memo8: TMemo;

Memo9: TMemo;

Memo10: TMemo;

StaticText11: TStaticText;

StaticText12: TStaticText;

StaticText13: TStaticText;

procedure Button1Click (Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

Form1: TForm1;

implementation

uses Okno2;

{$R *.dfm}

procedure TForm1. Button1Click (Sender: TObject);

const

h=0.001;

eps=0.01;

alfa=9;

k=10;

Pk=0.1;

kk=1300;

t=1.3;

var

Nk:array [1.10] of longint;

Dk:array [1.10] of real;

ntreb, n: real;

i, j, jj, jj1: integer;

f:boolean;

Dx, Mx, Mxcp, Sx, Sx2, x, qk2, qk, g, a, Disp, treal, xn: real;

begin

randomize;

ntreb:=200;

n:=0;

f:=true;

jj1:=0;

Memo1.Lines.Clear;

Memo2.Lines.Clear;

Memo3.Lines.Clear;

Memo4.Lines.Clear;

Memo5.Lines.Clear;

Memo6.Lines.Clear;

Memo7.Lines.Clear;

Memo8.Lines.Clear;

Memo9.Lines.Clear;

Memo10.Lines.Clear;

Memo11.Lines.Clear;

Memo12.Lines.Clear;

while n

Dx:=0;

Mx:=0;

qk:=0;

Mxcp:=0;

Disp:=0;

for i:=1 to k do begin

Sx:=0;

Sx2:=0;

if f then Nk[i]: =9

else

Nk[i]:=round (ntreb*Pk*sqrt (Dk[i])/qk2);

for j:=1 to Nk[i] do begin

g:=random*0.4+1;

a:=random*Pk*0.2+0.6+Pk*(i-1)*0.2;

x:=1;

for jj:=1 to kk do begin x:=(a*x+g*a*1.2)*h+x;

end;

Sx:=Sx+x;

Sx2:=Sx2+sqr (x);

end;

Mx:=Sx/Nk[i];

Dk[i]:=Sx2/Nk[i]-sqr (Sx/Nk[i]);

Dx:=Dx+(Pk*Pk)*Dk[i]/Nk[i];

Mxcp:=Mxcp+Pk*Sx/Nk[i];

case i of

1: begin //Memo1.Lines.Clear;

Memo1.Lines.Add ('Mx='+FloatToStr (Mx));

Memo1.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo1.Lines.Add ('Disp='+FloatToStr (Disp));

Memo1.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo1.Lines.Add ('Sx='+FloatToStr (Sx));

Memo1.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo1.Lines.Add ('n='+FloatToStr (n));

Memo1.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo1.Lines.Add ('***********');

end;

2: begin //Memo2.Lines.Clear;

Memo2.Lines.Add ('Mx='+FloatToStr (Mx));

Memo2.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo2.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo2.Lines.Add ('Sx='+FloatToStr (Sx));

Memo2.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo2.Lines.Add ('n='+FloatToStr (n));

Memo2.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo2.Lines.Add ('***********');end;

3: begin //Memo3.Lines.Clear;

Memo3.Lines.Add ('Mx='+FloatToStr (Mx));

Memo3.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo3.Lines.Add ('Disp='+FloatToStr (Disp));

Memo3.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo3.Lines.Add ('Sx='+FloatToStr (Sx));

Memo3.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo3.Lines.Add ('n='+FloatToStr (n));

Memo3.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo3.Lines.Add ('***********');end;

4: begin// Memo4.Lines.Clear;

Memo4.Lines.Add ('Mx='+FloatToStr (Mx));

Memo4.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo4.Lines.Add ('n='+FloatToStr (n));

Memo4.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo4.Lines.Add ('Disp='+FloatToStr (Disp));

Memo4.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo4.Lines.Add ('Sx='+FloatToStr (Sx));

Memo4.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo4.Lines.Add ('***********');end;

5: begin //Memo5.Lines.Clear;

Memo5.Lines.Add ('Mx='+FloatToStr (Mx));

Memo5.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo5.Lines.Add ('n='+FloatToStr (n));

Memo5.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo5.Lines.Add ('Disp='+FloatToStr (Disp));

Memo5.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo5.Lines.Add ('Sx='+FloatToStr (Sx));

Memo5.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo5.Lines.Add ('***********');end;

6: begin //Memo6.Lines.Clear;

Memo6.Lines.Add ('Mx='+FloatToStr (Mx));

Memo6.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo6.Lines.Add ('n='+FloatToStr (n));

Memo6.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo6.Lines.Add ('Disp='+FloatToStr (Disp));

Memo6.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo6.Lines.Add ('Sx='+FloatToStr (Sx));

Memo6.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo6.Lines.Add ('***********');end;

7: begin //Memo7.Lines.Clear;

Memo7.Lines.Add ('Mx='+FloatToStr (Mx));

Memo7.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo7.Lines.Add ('n='+FloatToStr (n));

Memo7.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo7.Lines.Add ('Disp='+FloatToStr (Disp));

Memo7.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo7.Lines.Add ('Sx='+FloatToStr (Sx));

Memo7.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo7.Lines.Add ('***********');end;

8: begin //Memo8.Lines.Clear;

Memo8.Lines.Add ('Mx='+FloatToStr (Mx));

Memo8.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo8.Lines.Add ('n='+FloatToStr (n));

Memo8.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo8.Lines.Add ('Disp='+FloatToStr (Disp));

Memo8.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo8.Lines.Add ('Sx='+FloatToStr (Sx));

Memo8.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo8.Lines.Add ('***********');end;

9: begin //Memo9.Lines.Clear;

Memo9.Lines.Add ('Mx='+FloatToStr (Mx));

Memo9.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo9.Lines.Add ('n='+FloatToStr (n));

Memo9.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo9.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo9.Lines.Add ('Disp='+FloatToStr (Disp));

Memo9.Lines.Add ('Sx='+FloatToStr (Sx));

Memo9.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo9.Lines.Add ('***********');end;

10: begin //Memo10.Lines.Clear;

Memo10.Lines.Add ('Mx='+FloatToStr (Mx));

Memo10.Lines.Add ('Dx='+FloatToStr (Dk[i]));

Memo10.Lines.Add ('n='+FloatToStr (n));

Memo10.Lines.Add ('ntreb='+FloatToStr (ntreb));

Memo10.Lines.Add ('Disp='+FloatToStr (Disp));

Memo10.Lines.Add ('Nk='+FloatToStr (Nk[i]));

Memo10.Lines.Add ('Sx='+FloatToStr (Sx));

Memo10.Lines.Add ('Sx2='+FloatToStr (Sx2));

Memo10.Lines.Add ('***********');end;

end;

// Form2. Visible:=true;

qk:=Pk*sqrt (Dk[i])+qk;

Disp:=Disp+Dk[i]*qk/(Pk*sqrt (Dk[i]));

end;

qk2:=qk;

Memo11.Lines.Add ('Mxcp='+floattostr (Mxcp));

Memo11.Lines.Add ('Dxcp='+floattostr (Dx));

Memo11.Lines.Add ('***********');

edit1.Text:=floattostr (Mxcp);

edit2.Text:=floattostr (Dx);

n:=ntreb;

ntreb:=round (alfa*Disp*sqr (Pk)/sqr (eps));

Memo12.Lines.Add ('n='+floattostr (n));

Memo12.Lines.Add ('ntreb='+floattostr (ntreb));

Memo12.Lines.Add ('***********');

edit3.text:=floattostr (ntreb);

f:=false;

end.