Алтаев А.А., Бундаев В.В. Компьютерное моделирование. Лабораторный практикум

61

числа наблюдений. Если имеется два фактора и для каждого

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

эксперимента потребуется

S = 3

2

= 9 наблюдений, для

четырех факторов и четырех уровней -

S = 4

4

= 256

наблюдений. Для уменьшения числа наблюдений

применяют план эксперимента с варьированием всех

m

факторов на двух уровнях:

x

imin

, x

imax

. Количество

наблюдений функции отклика

в этом случае составит

S = 2

m

, а саму функцию

ψ

приближенно можно представить

в виде линейного полинома:

∑

≤≤≤≤≤≤≤≤≤≤

++++=

mkji

kjikij

mji

jiij

mi

ii

xxxbxxbxbby

...1

...

11

0

...Κ

(6.1)

Коэффициент b

i

отражает влияние фактора x

i

на

реакцию

y, коэффициент b

ij

– взаимовлияние факторов x

i

и

x

j

, коэффициент b

ij…k

– взаимовлияние факторов x

i

, x

j

, …,

x

k

.

В таблице 6.1 приведен план проведения

двухфакторного эксперимента:

S = 2

2

.

Таблица 6.1

План ПФЭ

Номер

испытания

x

1

x

2

Реакция

y

1

x

1min

x

2min

y

1

2

x

1min

x

2max

y

2

3

x

1max

x

2min

y

3

4

x

1max

x

2max

y

4

Функция отклика

ψ

при этом будет иметь вид

y = b

0

+ b

1

x

1

+ b

2

x

2

+ b

3

x

1

x

2

.

Для ПФЭ типа 2

3

план эксперимента приведен в

таблице 6.2.

62

Таблица 6.2.

План ПФЭ

Номер

испытания

x

1

x

3

x

3

Реакция

y

1

x

1min

x

2min

x

3min

y

1

2

x

1min

x

2min

x

3max

y

2

3

x

1min

x

2max

x

3min

y

3

4

x

1min

x

2max

x

3max

y

4

5

x

1max

x

2min

x

3min

y

5

6

x

1max

x

2min

x

3max

y

6

7

x

1max

x

2max

x

3min

y

7

8

x

1max

x

2max

x

3max

y

8

Уравнение регрессии

ψ

можно представить в виде

y = b

0

+ b

1

x

1

+ b

2

x

2

+ b

3

x

3

+ b

4

x

1

x

2

+ b

5

x

2

x

3

+ b

6

x

1

x

3

+

+ b

7

x

1

x

2

x

3

.

Коэффициент b

7

учитывает взаимовлияние всех трех

факторов друг на друга.

Основная идея регрессионного анализа

Для исследования взаимосвязи между величинами Х

(вход) и

y (выход) используются методы корреляционного и

регрессионного анализа. Результаты корреляционного

анализа позволяют сделать вывод о степени зависимости

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

методами регрессионного анализа.

На значение величины

y оказывают влияние

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

связи между величинами

Х и y определяется линией

63

регрессии, показывающей, как в

среднем изменяется

величина

y при изменении входной величины Х, т.е.

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

y

c X. Эту связь характеризуют условным математическим

ожиданием величины

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

величина

Х приняла определенное значение, а

аппроксимирующая функция строится

как функция

регрессии

ψ

(X, B)

≈

M[y/X], где B - неизвестные

параметры уравнения регрессии.

Задача регрессионного анализа ставится следующим

образом. Для каждого i-го опыта имеется набор значений

(

x

i1

,...,x

im

) входных параметров

x

1

÷ x

m

и соответствующее им значение выходного

параметра

y

i

. Пример опытных данных приведен в

таблице 6.3.

Таблица 6.3

Факторы

Номер

испытания

x

1

x

2

--- x

m

Измеренная в ходе

эксперимента

реакция

y

1

x

11

x

12

--- x

1m

y

1

2

x

21

x

22

--- x

2m

y

2

-----

----- ----- --- ----- -----

k x

k1

x

k2

--- x

km

y

k

Задача сводится к определению значений

коэффициентов уравнения регрессии

b

0

, b

1

, ... , b

n

, которые с

определенной степенью вероятности будут отражать

влияние аргументов

x

i1

,..., x

im

на y.

Для определения

b

k

используется метод

наименьших квадратов

(МНК), смысл которого сводится

к минимизации функции:

()

min

~

1

2

→−=

∑

=

k

i

ii

yyZ ,

(6.2)

64

где

y

i

- фактическое (экспериментальное) значение

выходной переменной,

~

y

i

- рассчитанное по уравнению регрессии значение

выходной переменной.

Выражение (6.2) с учетом (6.1) будет иметь вид:

()()

∑

=

→++++−=

k

i

imiiniii

xxxbxbxbbyZ

1

2

2122110

min......

.

(6.3)

В данной работе для вычисления

b

0

, b

1

,…, b

n

следует

использовать метод Нелдера-Мида, который реализован в

программе nm.pas.

Отметим, что современные системы моделирования

позволяют полностью автоматизировать весь процесс

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

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

проведения экспериментов и завершая оценкой влияния

факторов на функцию отклика. В качестве примера можно

привести общецелевую систему моделирования GPSS World

[3]. Исследователь для своей модели указывает функцию

отклика и факторы, назначает нижнюю и верхнюю границы

интервала варьирования каждого фактора, выбирает тип

эксперимента: дисперсионный или регрессионный анализ.

Дисперсионный анализ (отсеивающий эксперимент)

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

величину (отклик). Регрессионный анализ

(оптимизирующий эксперимент) позволяет получить

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

уравнения координаты точки экстремума.

65

Задание

Цель работы: определить коэффициенты b

0

, b

1

, ... ,

b

n

уравнения регрессии y=

ψ

(X) на основе обработки

результатов эксперимента.

В соответствии с вариантом составить GPSS-модель.

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

следующим указаниям:

− для каждого фактора определить по два уровня -

x

imin

, x

imax

;

− в Excel построить таблицу, аналогичную

таблице 6.3.

Провести эксперименты над GPSS-моделью и

полученные результаты

y

i

занести в таблицу 6.3.

Для вычисления коэффициентов

b

i

уравнения

регрессии формулу (6.3) внести в программу nm.pas. Для

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

как типизированную константу. Например,

const

ur=2; {

число уровней}

m=2;

{число факторов}

fact:array[1.. ur, 1..m] of real = …

Методика вычислений такая же как и в предыдущей

работе. Результаты расчетов занести в таблицу,

аналогичную 5.1.

После определения коэффициентов

b

i

следует

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

функции отклика, добавив к таблице 6.3. столбец

y

~

с

вычисленными по (6.1) теоретическими значениями.

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

эксперимент, построить график поверхности отклика. Для

этого заполнить таблицу 6.4.

66

Таблица 6.4

x

1min

x

1max

x

2min

y

1

y

3

x

2max

y

2

y

4

Тип графика выбрать

Поверхность, проволочная.

На основе анализа полученного регрессионного

уравнения сделайте вывод о том, какие меры следует

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

Варианты

Варьируемые факторы

Вариант

Задача

Функция отклика

Первый

фактор

Второй

фактор

1

коэффициент

использования первого

процессора

2

коэффициент

использования второго

процессора

3

коэффициент

использования внешней

памяти

4

среднее время нахождения

транзакта в системе

5,

6,

7

средняя длина очереди:

5) A1;

6) A2;

7) A3

8,

9,

10

1

максимальная длина

очереди:

8) A1;

9) A2;

10) A3

интервал моделирования

x

1min

=100 задач

x

1max

=500 задач

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

x

2min

=2 ед.вр.

x

2max

=8 ед.вр.

67

Варьируемые факторы

Вариант

Задача

Функция отклика

Первый

фактор

Второй

фактор

11

коэффициент

использования первого

процессора

12

коэффициент

использования второго

процессора

13

коэффициент

использования внешней

памяти

14

среднее время нахождения

транзакта в системе

15,

16,

17

средняя длина очереди:

15) A1;

16) A2;

17) A3

18,

19,

20

максимальная длина

очереди:

18) A1;

19) A2;

20) A3

21

В задаче 1 предусмотреть уход транзактов из системы без

обслуживания, если длина очереди А1 превысит 5

количество заявок,

покинувших систему без

обработки.

интервал моделирования

x

1min

=100 задач

x

1max

=500 задач

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

x

2min

=2 ед.вр.

x

2max

=8 ед.вр.

22,

23,

24

коэффициент

использования:

22) A1;

23) A2;

24) A3

25,

26

2

средняя длина очереди:

25) A1;

26) A4

единица времени

x

1min

= 1 мс

x

1max

= 1,7 мс

количество дисководов

x

2min

= 3

x

2max

= 4

68

Задача 1

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

двухпроцессорной ЭВМ с общей памятью, разделенной на

восемь блоков. Каждой задаче отводится при ее решении

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

распределены по экспоненциальному закону со средним

временем 8 единиц времени, время обработки порции

информации также подчинено экспоненциальному закону с

интенсивностью v1=5 в процессоре CPU1 и с v2=2 в

процессоре CPU2.

Между обработкой порций с вероятностью 0,6

возможно обращение к внешней памяти, в которой время

обслуживания распределено равномерно в диапазоне [2,8].

С вероятностью 0,4 задачи оказываются решенными и

покидают систему. Моделирование выполнить на отрезке

времени, соответствующем решению не менее 100 задач.

Ниже представлена схема имитационной модели и текст

программы на языке GPSS.

MEM STORAGE 8

EXP FUNCTION RN1,C12

0,0/.2,.22/.4,.51/.5,.69/.6,.92/.7,1.2/.8,1

.61/

.9,2.3/.95,3/.99,4.6/.999,6.9/1,100

GENERATE 8,FN$EXP,,100

QUEUE A1

ENTER MEM,1

DEPART A1

M6 QUEUE A2

TRANSFER BOTH,M1,M2

M1 SEIZE CPU1

DEPART A2

ADVANCE 5,FN$EXP

69

RELEASE CPU1

TRANSFER ,M3

M2 SEIZE CPU2

DEPART A2

ADVANCE 2,FN$EXP

RELEASE CPU2

M3 TRANSFER .6,M5,M4

M4 QUEUE A3

SEIZE DISK

DEPART A3

ADVANCE 5,3

RELEASE DISK

TRANSFER ,M6

M5 LEAVE MEM,1

TERMINATE 1

Задача 2

Промоделировать работу устройства дисковой

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

Запросы поступают равновероятные ко всем дисководам.

Обработка запроса включает установку головки (при этом

канал не требуется) и обмен данными через канал.

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

распределены по экспоненциальному закону с v=6. Время

установки головки равномерно распределено в интервале 0 -

50 мс. Время обмена данными равно 1,7 мс (за единицу

времени принять 1,7 мс).

Ниже приведен текст программы на языке GPSS.

EXP FUNCTION RN1,C12

0,0/.2,.22/.4,.51/.5,.69/.6,.92/.7,1.2/.8,1

.61/

.9,2.3/.95,3/.99,4.6/.999,6.9/1,100

GENERATE 6,FN$EXP

70

TRANSFER .333,M2,M1

M2 TRANSFER .5,M4,M3

M1 QUEUE A1

SEIZE DISK1

DEPART A1

ASSIGN 1,DISK1

ADVANCE 15,15

TRANSFER ,M5

M3 QUEUE A2

SEIZE DISK2

DEPART A2

ASSIGN 1,DISK2

ADVANCE 15,15

TRANSFER ,M5

M4 QUEUE A3

SEIZE DISK3

DEPART A3

ASSIGN 1,DISK3

ADVANCE 15,15

M5 QUEUE A4

SEIZE CAN

DEPART A4

ADVANCE 1

RELEASE CAN

RELEASE P1

TERMINATE 1

71

Литература

1. Банди Б. Методы оптимизации. Вводный курс: Пер. с

англ. – М.: Радио и связь, 1988.

2.

Боев В.Д. Моделирование систем. Инструментальные

средства GPSS World: Учеб. пособие. – СПб.:

БХВ-Петербург, 2004. – 368 с.

3.

Голованов О. В., Дуваков С., Смирнов В. Н.

Моделирование сложных дискретных систем на ЭВМ

третьего поколения (опыт применения GPSS). – М.:

Энергия, 1978. – 160 с.

4.

Курс лекций по дисциплине «Компьютерное

моделирование»/ Сост. Алтаев А.А. – Улан-Удэ, Изд-во

ВСГТУ, 2001. – 63 с.

5.

Наставление по GPSS/PC. Minuteman Software: Пер. с

англ. под ред. И. М. Якимова.– Казань, 1997. — 320 с.

6.

Советов Б. Я., Яковлев С. А. Моделирование систем:

Учебник для ВУЗов. – М.: Высшая школа, 1999.

7.

Советов Б. Я., Яковлев С. А. Моделирование систем.

Лабораторный практикум. – М.: Высшая школа, 1989.

8.

Федоров В.Н. Моделирование дискретных систем.

Учебное пособие – М.: МГАПИ, 2005. – 92 с.

9.

Шеннон Р. Дж. Имитационное моделирование систем -

искусство и наука. – М.: Мир, 1978. - 418 с.

10.

Шрайбер Т. Дж. Моделирование на GPSS. – М.:

Машиностроение, 1980. - 592 с.

72

Приложение

Метод деформируемого многогранника

Метод деформируемого многогранника или, второе

его название, - метод Нелдера-Мида (Nelder-Mead)

относится к методам оптимизации нулевого порядка, то

есть к методам, использующим в ходе поиска минимума

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

для поиска минимума помимо значений самой целевой

функции используют еще и значения ее первой

производной, методы второго порядка – значения второй

производной и т.д.

Рассмотрим алгоритм метода в его классическом

варианте. Итак, пусть необходимо найти минимум целевой

функции от

n переменных. В n-мерном эвклидовом

пространстве строится многогранник с

n+1 вершиной.

Например, в случае одной переменной многогранником

x

1

x

1

Рис. П.1

x

2

а

б

73

будет отрезок (рис. П.1,а), для двух переменных

−

треугольник (рис. П.1,б), при

n=3 − это тетраэдр.

Далее используются следующие обозначения:

x

ij

−

координата

i-й вершины, j − номер переменной, f

i

—

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

i-й вершины. Помимо матрицы

X(N+1,N) и вектора F(N+1) применяются вспомогательные

векторы

XO(N), XC(N), XR(N), XE(N).

В качестве исходных данных задаются:

− n − число переменных;

− координаты первой вершины многогранника

− x

1

j, где nj ,1= ;

− начальный размер многогранника ∆.

Координаты других вершин определяются по

формуле







+=∆+

+≠

=

1,

1

jiеслиx

x

j

ij

, 1,2 += ni .

При этом создается равносторонний многогранник,

называемый симплексом. Сам метод по этой причине имеет

третье название

− симплексный метод (не путать с

симплекс-методом, который используется при решении

задач линейного программирования).

Далее вычисляются значения целевой функции в

вершинах

− f

i

.

Дальнейшие действия состоят в перемещении

построенного многогранника в область минимума функции.

Алгоритм перемещения рассмотрим на примере целевой

функции от двух переменных (рис. П.2).

1. Определяется вершина с наибольшим значением

целевой функции. Номер такой вершины обозначим через

H. Таким образом, f

h

= max{f

1

, f

2

,…, f

n

}. Также определяется

вершина

L с наименьшей функцией: f

l

= min{f

1

, f

2

,…, f

n

}.

Комментарий по алгоритму. Основная идея метода

состоит в том, что на каждом шаге очередная «наихудшая»

74

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

целевой функции будет меньше. Это позволит смещать

многогранник в направлении минимума.

2. Вычисляются координаты центра тяжести

O всех

вершин многогранника, исключая

H :

n

xx

xo

hj

n

i

ij

j

−













=

∑

+1

.

3. Вычисляются координаты точки

R, лежащей на

продолжении отрезка

HO, причем длины отрезков HO и OR

равны:

xr

j

=2·xo

j

−

x

hj

.

4. Вычисляется целевая функция в точке R − f

r

.

5. Если

f

r

< f

h

, то идти к п. 6, иначе к п. 13.

x

1

R

E

O

C

H

Рис. П.2

x

2

Линии

равного

уровня

функции

75

Комментарий по п.5. Если выполняется условие f

r

< f

h

, то это означает, что точка R вполне может быть

принята в качестве новой вершины. Операция замены

H на

R называется отражением. Но имеет смысл проверить и

более отдаленную точку

− E, выбор которой позволит

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

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

многогранника, что повышает эффективность алгоритма в

случае, если размеры исходного многогранника были

выбраны небольшими и он находится далеко от точки

минимума.

6. Проверяется точка

E в качестве альтернативы

точке

R:

xe

j

=2·x

rj

−

xo

j

.

7. Вычисляется f

e

.

8. Если

f

e

< f

h

, то идти к п. 9, иначе к п. 11.

9. В качестве новой вершины выбирается точка

E

(выполняется

растяжение):

x

hj

=xe

j

.

10. Идти к п. 20.

11. В качестве новой вершины выбирается точка

R

(выполняется

отражение):

x

hj

=xr

j

.

12. Идти к п. 20.

13. Так как точка

R не подходит в качестве новой

вершины (

f

r

> f

h

), то следующей рассматривается точка C,

лежащая посередине отрезка

HO:

2

jhj

j

xox

xc

+

=

.

14. Вычисляется

f

c.

15. Если f

c

< f

h

, то идти к п. 16, иначе к п. 18.

16. В качестве новой вершины выбирается точка

C

(выполняется

сжатие):

76

x

hj

=xc

j

.

17. Идти к п. 20.

18. Так как точки

R и C не подошли в качестве

альтернативы вершине

H, выполняется всеобщее сжатие

многогранника, причем наилучшая вершина

L остается на

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

они сместились в сторону наилучшей вершины:

.,1,1,1,

2

njni

xx

x

ljij

ij

=+=

+

=

19. Вычисляются функции

f

i

.

20. Выводятся на печать:

l − номер наилучшей

вершины,

f

l

− значение целевой функции в данной вершине,

x

lj

− координаты вершины.

21. Проверяется условие выхода из цикла. Для этого

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

функций для всех вершин, которое сравнивается с

погрешностью расчета

ε.

1

2

1

2

+













−

=

∑

+

=

+

=

n

f

n

i

n

i

σ

.

Комментарий по п.21. После перемещения

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

операции сжатия и вершины постепенно стягиваются в одну

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

значениям. Это обстоятельство используется для

завершения вычислений.

22. Если

σ < ε, то идти к п. 23, иначе к п. 1.

23. Выводятся на печать:

f

l

− значение целевой

функции в наилучшей вершине,

x

lj

− координаты вершины.

77

Комментарий по алгоритму. Поскольку алгоритм

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

и формы многогранника за счет растяжения и сжатия, метод

получил название

деформируемого многогранника.

Ниже дана распечатка текста программы, в которой

реализован алгоритм несколько отличающийся от

рассмотренного выше. В частности, для построения

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

чисел. Имеются и другие отличия.

program nm;

const

nMax=10;

type dim= array[1..nMax+1,1..nMax] of real;

vect1=array [1..nMax] of real;

vect2=array [1..nMax+1] of real;

var

s: dim;

x,xh,xg,xl,xo,xr,xc,xe:vect1;

f: vect2;

ko,z,k,fl,fh,s1,s2,sig:real;

j,n,jc,i,tev,l,u,u1,h,g:integer;

fg,fr,fc,fe:real;

{Выбор вершины}

procedure choice(fn:real; x:vect1);

var j: integer;

begin

for j:=1 to n do

s[h,j]:=x[j];

f[h]:=fn

end;

{Вычисление целевой функции}

procedure functmin(var z:real; x:vect1);

begin

tev:=tev+1;

z:=100*sqr(x[2]-x[1])+sqr(1-x[1]);

78

end;

begin

randomize;

writeln('симплексный метод');

writeln('введи число пеpеменных');

readln(n);

tev:=0;

ko:=1;

writeln('Начальное пpиближение (x><0)');

for j:=1 to n do

begin

writeln('введите x',j);

read(s[1,j]);

if s[1,j]=0 then

s[1,j]:=0.001

end;

writeln('введите длину шага ( >1 )');

readln(k);

if k=1 then k:=1.1;

u1:=n+trunc(1.2*ln(k));

fl:=1e+20;

l:=1;

for i:=1 to n+1 do

begin

f[i]:=1e+20;

for u:=1 to u1 do

begin

for j:=1 to n do

x[j]:=s[l,j]*exp((1-

2*random)*ln(k));

functmin(z,x);

ko:=0.5*k+0.5;

if f[i]>z then

begin

f[i]:=z;

for j:=1 to n do

s[i,j]:=x[j];

if z<fl then

79

begin

l:=i;

fl:=f[i];

ko:=k*3-2;

k:=ko;

end;

repeat

fh:=-1e+20;

fl:=1e+20;

for i:=1 to n+1 do

begin

if f[i]>fh then

begin

fh:=f[i];

h:=i;

end;

if f[i]<fl then

begin

fl:=f[i];

l:=i;

end;

fg:=-1e+20;

for i:=1 to n+1 do

begin

if i<>h then

if f[i]>fg then

begin

fg:=f[i];

g:=i;

end;

for j:=1 to n do

begin

xo[j]:=0;

80

for i:=1 to n+1 do

xo[j]:=xo[j]+s[i,j];

xo[j]:=(xo[j]-s[h,j])/n;

xh[j]:=s[h,j];

xg[j]:=s[g,j];

xl[j]:=s[l,j];

end;

for j:=1 to n do

xr[j]:=2*xo[j]-xh[j];

functmin(fr,xr);

if fr<fl then

begin

for j:=1 to n do

xe[j]:=2*xr[j]-xo[j];

functmin(fe,xe);

if fe<fr then choice(fe,xe)

else choice(fr,xr)

end

else

begin

if fr>fg then

begin

if fr<fh then choice(fr,xr)

else

begin

for j:=1 to n do

xc[j]:=(xh[j]+xo[j])/2;

functmin(fc,xc);

if fc>fh then

begin {}

for i:=1 to n+1 do

begin

for j:=1 to n do

begin

s[i,j]:=(s[i,j]+xl[j])/2;

x[j]:=s[i,j];

end;

functmin(f[i],x);

Алтаев А.А., Бундаев В.В. Компьютерное моделирование. Лабораторный практикум

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